SyFi 0.3
|
00001 // Copyright (C) 2006-2009 Kent-Andre Mardal and Simula Research Laboratory. 00002 // Licensed under the GNU GPL Version 2, or (at your option) any later version. 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 // remove previously computed basis functions and dofs 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 // FIXME check Andrew/Ola numbering 00215 equations.append( dofvdzy == GiNaC::numeric(0)); 00216 // FIXME check Andrew/Ola numbering 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 // matrix of symbols 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 // identity 00247 const GiNaC::ex _ex1(1); 00248 for (unsigned i=0; i<ncols; ++i) 00249 id(i, i) = _ex1; 00250 00251 // invert the matrix 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 } // namespace SyFi