SyFi 0.3
CrouzeixRaviart.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 "CrouzeixRaviart.h"
00005 
00006 #include "tools.h"
00007 
00008 using std::cout;
00009 using std::endl;
00010 using std::string;
00011 
00012 namespace SyFi
00013 {
00014 
00015         CrouzeixRaviart:: CrouzeixRaviart() : StandardFE()
00016         {
00017                 order = 1;
00018         }
00019 
00020         CrouzeixRaviart:: CrouzeixRaviart (Polygon& p, unsigned int order) : StandardFE(p, order)
00021         {
00022                 compute_basis_functions();
00023         }
00024 
00025         void CrouzeixRaviart:: compute_basis_functions()
00026         {
00027 
00028                 // remove previously computed basis functions and dofs
00029                 Ns.clear();
00030                 dofs.clear();
00031 
00032                 if ( p == NULL )
00033                 {
00034                         throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00035                 }
00036 
00037                 if (order != 1)
00038                 {
00039                         throw(std::logic_error("Only Crouziex-Raviart elements of order 1 is possible"));
00040                 }
00041 
00042                 // see e.g. Brezzi and Fortin book page 116 for the definition
00043 
00044                 if ( p->str().find("ReferenceLine") != string::npos )
00045                 {
00046                         cout <<"Can not define the Raviart-Thomas element on a line"<<endl;
00047                 }
00048                 else if ( p->str().find("Triangle") != string::npos )
00049                 {
00050 
00051                         description = "CrouzeixRaviart_2D";
00052 
00053                         Triangle& triangle = (Triangle&)(*p);
00054 
00055                         // create the polynomial space
00056                         GiNaC::ex polynom_space = bernstein(1, triangle, "a");
00057                         GiNaC::ex polynom = polynom_space.op(0);
00058                         GiNaC::lst variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00059                         GiNaC::ex basis = polynom_space.op(2);
00060 
00061                         // create the dofs
00062                         GiNaC::symbol t("t");
00063                         for (int j=0; j< 3; j++)
00064                         {
00065 
00066                                 // solve the linear system to compute
00067                                 // each of the basis functions
00068                                 GiNaC::lst equations;
00069                                 for (int i=0; i< 3; i++)
00070                                 {
00071 
00072                                         Line line = triangle.line(i);
00073                                         //                          GiNaC::ex dofi = line.integrate(polynom);
00074                                         GiNaC::lst midpoint = GiNaC::lst(
00075                                                 (line.vertex(0).op(0) +  line.vertex(1).op(0))/2,
00076                                                 (line.vertex(0).op(1) +  line.vertex(1).op(1))/2);
00077                                         dofs.insert(dofs.end(), midpoint);
00078 
00079                                         //                          GiNaC::ex dofi = polynom.subs( x == midpoint.op(0)).subs( y == midpoint.op(1));
00080                                         GiNaC::ex dofi = line.integrate(polynom);
00081                                         equations.append( dofi == dirac(i,j));
00082 
00083                                         if (j == 1)
00084                                         {
00085                                                 //                                  GiNaC::lst d = GiNaC::lst(line.vertex(0) ,  line.vertex(1));
00086                                                 dofs.insert(dofs.end(), midpoint);
00087                                         }
00088 
00089                                 }
00090                                 GiNaC::ex sub = lsolve(equations, variables);
00091                                 GiNaC::ex Ni = polynom.subs(sub);
00092                                 Ns.insert(Ns.end(),Ni);
00093                         }
00094 
00095                 }
00096                 else if ( p->str().find("Tetrahedron") != string::npos )
00097                 {
00098 
00099                         description = "CrouzeixRaviart_3D";
00100 
00101                         Tetrahedron& tetrahedron = (Tetrahedron&)(*p);
00102                         GiNaC::ex polynom_space = bernstein(1, tetrahedron, "a");
00103                         GiNaC::ex polynom = polynom_space.op(0);
00104                         GiNaC::lst variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00105                         GiNaC::ex basis = polynom_space.op(2);
00106 
00107                         GiNaC::ex bernstein_pol;
00108 
00109                         GiNaC::symbol t("t");
00110                         // dofs related to edges
00111                         for (int j=0; j< 4; j++)
00112                         {
00113 
00114                                 GiNaC::lst equations;
00115                                 for (int i=0; i< 4; i++)
00116                                 {
00117                                         Triangle triangle = tetrahedron.triangle(i);
00118                                         GiNaC::lst midpoint = GiNaC::lst(
00119                                                 (triangle.vertex(0).op(0) + triangle.vertex(1).op(0) + triangle.vertex(2).op(0))/3,
00120                                                 (triangle.vertex(0).op(1) + triangle.vertex(1).op(1) + triangle.vertex(2).op(1))/3,
00121                                                 (triangle.vertex(0).op(2) + triangle.vertex(1).op(2) + triangle.vertex(2).op(2))/3
00122                                                 );
00123 
00124                                         GiNaC::ex dofi = polynom.subs(x == midpoint.op(0)).subs(y == midpoint.op(1)).subs(z == midpoint.op(2));
00125                                         equations.append( dofi == dirac(i,j));
00126 
00127                                         if ( j == 1 )
00128                                         {
00129                                                 //                                  GiNaC::lst d = GiNaC::lst(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2));
00130                                                 dofs.insert(dofs.end(), midpoint);
00131                                         }
00132                                 }
00133                                 GiNaC::ex sub = lsolve(equations, variables);
00134                                 GiNaC::ex Ni = polynom.subs(sub);
00135                                 Ns.insert(Ns.end(),Ni);
00136 
00137                         }
00138                 }
00139         }
00140 
00141         // ------------VectorCrouzeixRaviart ---
00142 
00143         VectorCrouzeixRaviart:: VectorCrouzeixRaviart (Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order)
00144         {
00145                 size = size_ < 0 ? nsd: size_;
00146                 compute_basis_functions();
00147         }
00148 
00149         VectorCrouzeixRaviart:: VectorCrouzeixRaviart() : StandardFE()
00150         {
00151                 description = "VectorCrouzeixRaviart";
00152                 order = 1;
00153         }
00154 
00155         void VectorCrouzeixRaviart:: compute_basis_functions()
00156         {
00157 
00158                 if (order != 1)
00159                 {
00160                         throw(std::logic_error("Only Crouziex-Raviart elements of order 1 is possible"));
00161                 }
00162 
00163                 CrouzeixRaviart fe;
00164                 fe.set_polygon(*p);
00165                 fe.compute_basis_functions();
00166 
00167                 description = "Vector" + fe.str();
00168 
00169                 GiNaC::lst zero_list;
00170                 for (unsigned int s=1; s<= size ; s++)
00171                 {
00172                         zero_list.append(0);
00173                 }
00174 
00175                 for (unsigned int s=0; s< size ; s++)
00176                 {
00177                         for (unsigned int i=0; i< fe.nbf() ; i++)
00178                         {
00179                                 GiNaC::lst Nis = zero_list;
00180                                 Nis.let_op(s) = fe.N(i);
00181                                 GiNaC::ex Nmat = GiNaC::matrix(size,1,Nis);
00182                                 Ns.insert(Ns.end(), Nmat);
00183 
00184                                 GiNaC::lst dof = GiNaC::lst(fe.dof(i), s) ;
00185                                 dofs.insert(dofs.end(), dof);
00186                         }
00187                 }
00188         }
00189 
00190         void VectorCrouzeixRaviart:: set_size(unsigned int size_)
00191         {
00192                 size = size_;
00193         }
00194 
00195 }                                                                // namespace SyFi
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines