|
SyFi
0.3
|
|
SyFi
0.3
|
00001 // Copyright (C) 2006-2009 Kent-Andre Mardal and Simula Research Laboratory 00002 // 00003 // This file is part of SyFi. 00004 // 00005 // SyFi is free software: you can redistribute it and/or modify 00006 // it under the terms of the GNU General Public License as published by 00007 // the Free Software Foundation, either version 2 of the License, or 00008 // (at your option) any later version. 00009 // 00010 // SyFi is distributed in the hope that it will be useful, 00011 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 // GNU General Public License for more details. 00014 // 00015 // You should have received a copy of the GNU General Public License 00016 // along with SyFi. If not, see <http://www.gnu.org/licenses/>. 00017 00018 #include "Hermite.h" 00019 #include "tools.h" 00020 00021 using std::cout; 00022 using std::endl; 00023 using std::string; 00024 00025 namespace SyFi 00026 { 00027 00028 Hermite:: Hermite() : StandardFE() 00029 { 00030 description = "Hermite"; 00031 } 00032 00033 Hermite:: Hermite(Polygon& p, int order) : StandardFE(p,order) 00034 { 00035 compute_basis_functions(); 00036 } 00037 00038 void Hermite:: compute_basis_functions() 00039 { 00040 00041 // remove previously computed basis functions and dofs 00042 Ns.clear(); 00043 dofs.clear(); 00044 00045 if ( p == NULL ) 00046 { 00047 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00048 } 00049 00050 GiNaC::ex polynom_space; 00051 GiNaC::ex polynom; 00052 GiNaC::lst variables; 00053 GiNaC::lst equations; 00054 00055 if ( p->str().find("Line") != string::npos ) 00056 { 00057 00058 description = "Hermite_1D"; 00059 00060 polynom_space = legendre(3, 1, "a"); 00061 polynom = polynom_space.op(0); 00062 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00063 00064 for (int i=0; i< 2; i++) 00065 { 00066 GiNaC::ex v = p->vertex(i); 00067 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0))); 00068 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0))); 00069 equations.append( dofv == GiNaC::numeric(0)); 00070 equations.append( dofvdx == GiNaC::numeric(0)); 00071 00072 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), 0)); 00073 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), 1)); 00074 } 00075 00076 } 00077 00078 if ( p->str().find("Triangle") != string::npos ) 00079 { 00080 00081 description = "Hermite_2D"; 00082 00083 polynom_space = pol(3, 2, "a"); 00084 polynom = polynom_space.op(0); 00085 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00086 00087 for (int i=0; i<= 2; i++) 00088 { 00089 GiNaC::ex v = p->vertex(i); 00090 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0), y == v.op(1))); 00091 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0), y == v.op(1))); 00092 GiNaC::ex dofvdy = diff(polynom,y).subs(GiNaC::lst(x == v.op(0), y == v.op(1))); 00093 00094 equations.append( dofv == GiNaC::numeric(0)); 00095 equations.append( dofvdx == GiNaC::numeric(0)); 00096 equations.append( dofvdy == GiNaC::numeric(0)); 00097 00098 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 0)); 00099 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 1)); 00100 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 2)); 00101 } 00102 GiNaC::ex midpoint = GiNaC::lst((p->vertex(0).op(0) + p->vertex(1).op(0) + p->vertex(2).op(0))/3, 00103 (p->vertex(0).op(1) + p->vertex(1).op(1) + p->vertex(2).op(1))/3); 00104 GiNaC::ex dofm = polynom.subs(GiNaC::lst(x == midpoint.op(0), y == midpoint.op(1))); 00105 dofs.insert(dofs.end(), midpoint ); 00106 equations.append( dofm == GiNaC::numeric(0)); 00107 00108 } 00109 00110 else if ( p->str().find("Rectangle") != string::npos ) 00111 { 00112 00113 description = "Hermite_2D"; 00114 00115 polynom_space = legendre(3, 2, "a"); 00116 polynom = polynom_space.op(0); 00117 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00118 00119 for (int i=0; i< 4; i++) 00120 { 00121 GiNaC::ex v = p->vertex(i); 00122 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0), y == v.op(1))); 00123 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0), y == v.op(1))); 00124 GiNaC::ex dofvdy = diff(polynom,y).subs(GiNaC::lst(x == v.op(0), y == v.op(1))); 00125 GiNaC::ex dofvdyx = diff(diff(polynom,y),x).subs(GiNaC::lst(x == v.op(0), y == v.op(1))); 00126 equations.append( dofv == GiNaC::numeric(0)); 00127 equations.append( dofvdx == GiNaC::numeric(0)); 00128 equations.append( dofvdy == GiNaC::numeric(0)); 00129 equations.append( dofvdyx == GiNaC::numeric(0)); 00130 00131 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 0)); 00132 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 1)); 00133 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 2)); 00134 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 3)); 00135 } 00136 00137 } 00138 else if ( p->str().find("Tetrahedron") != string::npos ) 00139 { 00140 00141 description = "Hermite_3D"; 00142 00143 polynom_space = pol(3, 3, "a"); 00144 polynom = polynom_space.op(0); 00145 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00146 00147 for (int i=0; i<= 3; i++) 00148 { 00149 GiNaC::ex v = p->vertex(i); 00150 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00151 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00152 GiNaC::ex dofvdy = diff(polynom,y).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00153 GiNaC::ex dofvdz = diff(polynom,z).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00154 00155 equations.append( dofv == GiNaC::numeric(0)); 00156 equations.append( dofvdx == GiNaC::numeric(0)); 00157 equations.append( dofvdy == GiNaC::numeric(0)); 00158 equations.append( dofvdz == GiNaC::numeric(0)); 00159 00160 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 0)); 00161 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 1)); 00162 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 2)); 00163 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), 3)); 00164 00165 } 00166 GiNaC::ex midpoint1 = GiNaC::lst( 00167 (p->vertex(0).op(0)*2 + p->vertex(1).op(0) + p->vertex(2).op(0) + p->vertex(3).op(0))/5, 00168 (p->vertex(0).op(1)*2 + p->vertex(1).op(1) + p->vertex(2).op(1) + p->vertex(3).op(1))/5, 00169 (p->vertex(0).op(2)*2 + p->vertex(1).op(2) + p->vertex(2).op(2) + p->vertex(3).op(2))/5); 00170 00171 GiNaC::ex midpoint2 = GiNaC::lst( 00172 (p->vertex(0).op(0) + p->vertex(1).op(0)*2 + p->vertex(2).op(0) + p->vertex(3).op(0))/5, 00173 (p->vertex(0).op(1) + p->vertex(1).op(1)*2 + p->vertex(2).op(1) + p->vertex(3).op(1))/5, 00174 (p->vertex(0).op(2) + p->vertex(1).op(2)*2 + p->vertex(2).op(2) + p->vertex(3).op(2))/5); 00175 00176 GiNaC::ex midpoint3 = GiNaC::lst( 00177 (p->vertex(0).op(0) + p->vertex(1).op(0) + p->vertex(2).op(0)*2 + p->vertex(3).op(0))/5, 00178 (p->vertex(0).op(1) + p->vertex(1).op(1) + p->vertex(2).op(1)*2 + p->vertex(3).op(1))/5, 00179 (p->vertex(0).op(2) + p->vertex(1).op(2) + p->vertex(2).op(2)*2 + p->vertex(3).op(2))/5); 00180 00181 GiNaC::ex midpoint4 = GiNaC::lst( 00182 (p->vertex(0).op(0) + p->vertex(1).op(0) + p->vertex(2).op(0) + p->vertex(3).op(0)*2)/5, 00183 (p->vertex(0).op(1) + p->vertex(1).op(1) + p->vertex(2).op(1) + p->vertex(3).op(1)*2)/5, 00184 (p->vertex(0).op(2) + p->vertex(1).op(2) + p->vertex(2).op(2) + p->vertex(3).op(2)*2)/5); 00185 00186 GiNaC::ex dofm1 = polynom.subs(GiNaC::lst(x == midpoint1.op(0), y == midpoint1.op(1), z == midpoint1.op(2))); 00187 GiNaC::ex dofm2 = polynom.subs(GiNaC::lst(x == midpoint2.op(0), y == midpoint2.op(1), z == midpoint2.op(2))); 00188 GiNaC::ex dofm3 = polynom.subs(GiNaC::lst(x == midpoint3.op(0), y == midpoint3.op(1), z == midpoint3.op(2))); 00189 GiNaC::ex dofm4 = polynom.subs(GiNaC::lst(x == midpoint4.op(0), y == midpoint4.op(1), z == midpoint4.op(2))); 00190 00191 dofs.insert(dofs.end(), midpoint1 ); 00192 dofs.insert(dofs.end(), midpoint2 ); 00193 dofs.insert(dofs.end(), midpoint3 ); 00194 dofs.insert(dofs.end(), midpoint4 ); 00195 00196 equations.append( dofm1 == GiNaC::numeric(0)); 00197 equations.append( dofm2 == GiNaC::numeric(0)); 00198 equations.append( dofm3 == GiNaC::numeric(0)); 00199 equations.append( dofm4 == GiNaC::numeric(0)); 00200 00201 } 00202 else if ( p->str().find("Box") != string::npos ) 00203 { 00204 00205 description = "Hermite_3D"; 00206 00207 polynom_space = legendre(3, 3, "a"); 00208 polynom = polynom_space.op(0); 00209 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00210 00211 for (int i=0; i<= 7; i++) 00212 { 00213 GiNaC::ex v = p->vertex(i); 00214 GiNaC::ex dofv = polynom.subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00215 GiNaC::ex dofvdx = diff(polynom,x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00216 GiNaC::ex dofvdy = diff(polynom,y).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00217 GiNaC::ex dofvdyx = diff(diff(polynom,y),x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00218 GiNaC::ex dofvdz = diff(polynom,z).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00219 GiNaC::ex dofvdzy = diff(diff(polynom,z),y).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00220 GiNaC::ex dofvdzx = diff(diff(polynom,z),x).subs(GiNaC::lst(x == v.op(0), y == v.op(1), z == v.op(2) )); 00221 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) )); 00222 00223 equations.append( dofv == GiNaC::numeric(0)); 00224 equations.append( dofvdx == GiNaC::numeric(0)); 00225 equations.append( dofvdy == GiNaC::numeric(0)); 00226 equations.append( dofvdyx == GiNaC::numeric(0)); 00227 equations.append( dofvdz == GiNaC::numeric(0)); 00228 // FIXME check Andrew/Ola numbering 00229 equations.append( dofvdzy == GiNaC::numeric(0)); 00230 // FIXME check Andrew/Ola numbering 00231 equations.append( dofvdzx == GiNaC::numeric(0)); 00232 equations.append( dofvdzyx == GiNaC::numeric(0)); 00233 00234 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 0)); 00235 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 1)); 00236 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 2)); 00237 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 3)); 00238 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 4)); 00239 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 5)); 00240 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 6)); 00241 dofs.insert(dofs.end(), GiNaC::lst(v.op(0), v.op(1), v.op(2), 7)); 00242 00243 } 00244 00245 } 00246 00247 GiNaC::matrix b; GiNaC::matrix A; 00248 matrix_from_equations(equations, variables, A, b); 00249 00250 unsigned int ncols = A.cols(); 00251 GiNaC::matrix vars_sq(ncols, ncols); 00252 00253 // matrix of symbols 00254 for (unsigned r=0; r<ncols; ++r) 00255 for (unsigned c=0; c<ncols; ++c) 00256 vars_sq(r, c) = GiNaC::symbol(); 00257 00258 GiNaC::matrix id(ncols, ncols); 00259 00260 // identity 00261 const GiNaC::ex _ex1(1); 00262 for (unsigned i=0; i<ncols; ++i) 00263 id(i, i) = _ex1; 00264 00265 // invert the matrix 00266 GiNaC::matrix m_inv(ncols, ncols); 00267 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss); 00268 00269 for (unsigned int i=0; i<dofs.size(); i++) 00270 { 00271 b.let_op(i) = GiNaC::numeric(1); 00272 GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b)); 00273 00274 GiNaC::lst subs; 00275 for (unsigned int ii=0; ii<xx.nops(); ii++) 00276 { 00277 subs.append(variables.op(ii) == xx.op(ii)); 00278 } 00279 GiNaC::ex Nj= polynom.subs(subs).expand(); 00280 Ns.insert(Ns.end(), Nj); 00281 b.let_op(i) = GiNaC::numeric(0); 00282 } 00283 00284 } 00285 00286 } // namespace SyFi
1.7.6.1 
