SyFi 0.3
Hermite.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines