00001
00002
00003
00004 #include "Hermite.h"
00005 #include "tools.h"
00006
00007 using std::cout;
00008 using std::endl;
00009 using std::string;
00010
00011 namespace SyFi
00012 {
00013
00014 Hermite:: Hermite() : StandardFE()
00015 {
00016 description = "Hermite";
00017 }
00018
00019 Hermite:: Hermite(Polygon& p, int order) : StandardFE(p,order)
00020 {
00021 compute_basis_functions();
00022 }
00023
00024 void Hermite:: compute_basis_functions()
00025 {
00026
00027
00028 Ns.clear();
00029 dofs.clear();
00030
00031 if ( p == NULL )
00032 {
00033 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00034 }
00035
00036 GiNaC::ex polynom_space;
00037 GiNaC::ex polynom;
00038 GiNaC::lst variables;
00039 GiNaC::lst equations;
00040
00041 if ( p->str().find("Line") != string::npos )
00042 {
00043
00044 description = "Hermite_1D";
00045
00046 polynom_space = legendre(3, 1, "a");
00047 polynom = polynom_space.op(0);
00048 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00049
00050 for (int i=0; i< 2; i++)
00051 {
00052 GiNaC::ex v = p->vertex(i);
00053 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0)));
00054 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0)));
00055 equations.append( dofv == GiNaC::numeric(0));
00056 equations.append( dofvdx == GiNaC::numeric(0));
00057
00058 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), 0));
00059 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), 1));
00060 }
00061
00062 }
00063
00064 if ( p->str().find("Triangle") != string::npos )
00065 {
00066
00067 description = "Hermite_2D";
00068
00069 polynom_space = pol(3, 2, "a");
00070 polynom = polynom_space.op(0);
00071 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00072
00073 for (int i=0; i<= 2; i++)
00074 {
00075 GiNaC::ex v = p->vertex(i);
00076 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0), y == v.op(1)));
00077 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0), y == v.op(1)));
00078 GiNaC::ex dofvdy = diff(polynom,y).subs(GiNaC::lst(x == v.op(0), y == v.op(1)));
00079
00080 equations.append( dofv == GiNaC::numeric(0));
00081 equations.append( dofvdx == GiNaC::numeric(0));
00082 equations.append( dofvdy == GiNaC::numeric(0));
00083
00084 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 0));
00085 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 1));
00086 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 2));
00087 }
00088 GiNaC::ex midpoint = GiNaC::lst((p->vertex(0).op(0) + p->vertex(1).op(0) + p->vertex(2).op(0))/3,
00089 (p->vertex(0).op(1) + p->vertex(1).op(1) + p->vertex(2).op(1))/3);
00090 GiNaC::ex dofm = polynom.subs(GiNaC::lst(x == midpoint.op(0), y == midpoint.op(1)));
00091 dofs.insert(dofs.end(), midpoint );
00092 equations.append( dofm == GiNaC::numeric(0));
00093
00094 }
00095
00096 else if ( p->str().find("Rectangle") != string::npos )
00097 {
00098
00099 description = "Hermite_2D";
00100
00101 polynom_space = legendre(3, 2, "a");
00102 polynom = polynom_space.op(0);
00103 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00104
00105 for (int i=0; i< 4; i++)
00106 {
00107 GiNaC::ex v = p->vertex(i);
00108 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0), y == v.op(1)));
00109 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0), y == v.op(1)));
00110 GiNaC::ex dofvdy = diff(polynom,y).subs(GiNaC::lst(x == v.op(0), y == v.op(1)));
00111 GiNaC::ex dofvdyx = diff(diff(polynom,y),x).subs(GiNaC::lst(x == v.op(0), y == v.op(1)));
00112 equations.append( dofv == GiNaC::numeric(0));
00113 equations.append( dofvdx == GiNaC::numeric(0));
00114 equations.append( dofvdy == GiNaC::numeric(0));
00115 equations.append( dofvdyx == GiNaC::numeric(0));
00116
00117 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 0));
00118 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 1));
00119 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 2));
00120 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 3));
00121 }
00122
00123 }
00124 else if ( p->str().find("Tetrahedron") != string::npos )
00125 {
00126
00127 description = "Hermite_3D";
00128
00129 polynom_space = pol(3, 3, "a");
00130 polynom = polynom_space.op(0);
00131 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00132
00133 for (int i=0; i<= 3; i++)
00134 {
00135 GiNaC::ex v = p->vertex(i);
00136 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00137 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00138 GiNaC::ex dofvdy = diff(polynom,y).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00139 GiNaC::ex dofvdz = diff(polynom,z).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00140
00141 equations.append( dofv == GiNaC::numeric(0));
00142 equations.append( dofvdx == GiNaC::numeric(0));
00143 equations.append( dofvdy == GiNaC::numeric(0));
00144 equations.append( dofvdz == GiNaC::numeric(0));
00145
00146 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 0));
00147 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 1));
00148 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 2));
00149 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 3));
00150
00151 }
00152 GiNaC::ex midpoint1 = GiNaC::lst(
00153 (p->vertex(0).op(0)*2 + p->vertex(1).op(0) + p->vertex(2).op(0) + p->vertex(3).op(0))/5,
00154 (p->vertex(0).op(1)*2 + p->vertex(1).op(1) + p->vertex(2).op(1) + p->vertex(3).op(1))/5,
00155 (p->vertex(0).op(2)*2 + p->vertex(1).op(2) + p->vertex(2).op(2) + p->vertex(3).op(2))/5);
00156
00157 GiNaC::ex midpoint2 = GiNaC::lst(
00158 (p->vertex(0).op(0) + p->vertex(1).op(0)*2 + p->vertex(2).op(0) + p->vertex(3).op(0))/5,
00159 (p->vertex(0).op(1) + p->vertex(1).op(1)*2 + p->vertex(2).op(1) + p->vertex(3).op(1))/5,
00160 (p->vertex(0).op(2) + p->vertex(1).op(2)*2 + p->vertex(2).op(2) + p->vertex(3).op(2))/5);
00161
00162 GiNaC::ex midpoint3 = GiNaC::lst(
00163 (p->vertex(0).op(0) + p->vertex(1).op(0) + p->vertex(2).op(0)*2 + p->vertex(3).op(0))/5,
00164 (p->vertex(0).op(1) + p->vertex(1).op(1) + p->vertex(2).op(1)*2 + p->vertex(3).op(1))/5,
00165 (p->vertex(0).op(2) + p->vertex(1).op(2) + p->vertex(2).op(2)*2 + p->vertex(3).op(2))/5);
00166
00167 GiNaC::ex midpoint4 = GiNaC::lst(
00168 (p->vertex(0).op(0) + p->vertex(1).op(0) + p->vertex(2).op(0) + p->vertex(3).op(0)*2)/5,
00169 (p->vertex(0).op(1) + p->vertex(1).op(1) + p->vertex(2).op(1) + p->vertex(3).op(1)*2)/5,
00170 (p->vertex(0).op(2) + p->vertex(1).op(2) + p->vertex(2).op(2) + p->vertex(3).op(2)*2)/5);
00171
00172 GiNaC::ex dofm1 = polynom.subs(GiNaC::lst(x == midpoint1.op(0), y == midpoint1.op(1), z == midpoint1.op(2)));
00173 GiNaC::ex dofm2 = polynom.subs(GiNaC::lst(x == midpoint2.op(0), y == midpoint2.op(1), z == midpoint2.op(2)));
00174 GiNaC::ex dofm3 = polynom.subs(GiNaC::lst(x == midpoint3.op(0), y == midpoint3.op(1), z == midpoint3.op(2)));
00175 GiNaC::ex dofm4 = polynom.subs(GiNaC::lst(x == midpoint4.op(0), y == midpoint4.op(1), z == midpoint4.op(2)));
00176
00177 dofs.insert(dofs.end(), midpoint1 );
00178 dofs.insert(dofs.end(), midpoint2 );
00179 dofs.insert(dofs.end(), midpoint3 );
00180 dofs.insert(dofs.end(), midpoint4 );
00181
00182 equations.append( dofm1 == GiNaC::numeric(0));
00183 equations.append( dofm2 == GiNaC::numeric(0));
00184 equations.append( dofm3 == GiNaC::numeric(0));
00185 equations.append( dofm4 == GiNaC::numeric(0));
00186
00187 }
00188 else if ( p->str().find("Box") != string::npos )
00189 {
00190
00191 description = "Hermite_3D";
00192
00193 polynom_space = legendre(3, 3, "a");
00194 polynom = polynom_space.op(0);
00195 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00196
00197 for (int i=0; i<= 7; i++)
00198 {
00199 GiNaC::ex v = p->vertex(i);
00200 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00201 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00202 GiNaC::ex dofvdy = diff(polynom,y).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00203 GiNaC::ex dofvdyx = diff(diff(polynom,y),x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00204 GiNaC::ex dofvdz = diff(polynom,z).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00205 GiNaC::ex dofvdzy = diff(diff(polynom,z),y).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00206 GiNaC::ex dofvdzx = diff(diff(polynom,z),x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00207 GiNaC::ex dofvdzyx = diff(diff(diff(polynom,z),y),x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) ));
00208
00209 equations.append( dofv == GiNaC::numeric(0));
00210 equations.append( dofvdx == GiNaC::numeric(0));
00211 equations.append( dofvdy == GiNaC::numeric(0));
00212 equations.append( dofvdyx == GiNaC::numeric(0));
00213 equations.append( dofvdz == GiNaC::numeric(0));
00214
00215 equations.append( dofvdzy == GiNaC::numeric(0));
00216
00217 equations.append( dofvdzx == GiNaC::numeric(0));
00218 equations.append( dofvdzyx == GiNaC::numeric(0));
00219
00220 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 0));
00221 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 1));
00222 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 2));
00223 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 3));
00224 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 4));
00225 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 5));
00226 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 6));
00227 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 7));
00228
00229 }
00230
00231 }
00232
00233 GiNaC::matrix b; GiNaC::matrix A;
00234 matrix_from_equations(equations, variables, A, b);
00235
00236 unsigned int ncols = A.cols();
00237 GiNaC::matrix vars_sq(ncols, ncols);
00238
00239
00240 for (unsigned r=0; r<ncols; ++r)
00241 for (unsigned c=0; c<ncols; ++c)
00242 vars_sq(r, c) = GiNaC::symbol();
00243
00244 GiNaC::matrix id(ncols, ncols);
00245
00246
00247 const GiNaC::ex _ex1(1);
00248 for (unsigned i=0; i<ncols; ++i)
00249 id(i, i) = _ex1;
00250
00251
00252 GiNaC::matrix m_inv(ncols, ncols);
00253 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
00254
00255 for (unsigned int i=0; i<dofs.size(); i++)
00256 {
00257 b.let_op(i) = GiNaC::numeric(1);
00258 GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
00259
00260 GiNaC::lst subs;
00261 for (unsigned int ii=0; ii<xx.nops(); ii++)
00262 {
00263 subs.append(variables.op(ii) == xx.op(ii));
00264 }
00265 GiNaC::ex Nj= polynom.subs(subs).expand();
00266 Ns.insert(Ns.end(), Nj);
00267 b.let_op(i) = GiNaC::numeric(0);
00268 }
00269
00270 }
00271
00272 }