SyFi 0.3
BrezziDouglasMarini.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 "BrezziDouglasMarini.h"
00005 #include <fstream>
00006 
00007 #include "tools.h"
00008 
00009 using std::cout;
00010 using std::endl;
00011 using std::string;
00012 
00013 namespace SyFi
00014 {
00015 
00016         BrezziDouglasMarini:: BrezziDouglasMarini() : StandardFE()
00017         {
00018                 description = "BrezziDouglasMarini";
00019         }
00020 
00021         BrezziDouglasMarini:: BrezziDouglasMarini(Polygon& p, int order, bool pointwise_) : StandardFE(p, order)
00022         {
00023                 pointwise = pointwise_;
00024                 compute_basis_functions();
00025         }
00026 
00027         void BrezziDouglasMarini:: compute_basis_functions()
00028         {
00029 
00030                 if ( order < 1 )
00031                 {
00032                         throw(std::logic_error("Brezzi-Douglas-Marini elements must be of order 1 or higher."));
00033                 }
00034 
00035                 if ( p == NULL )
00036                 {
00037                         throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00038                 }
00039 
00040                 GiNaC::ex nsymb = GiNaC::symbol("n");
00041                 if (pointwise)
00042                 {
00043 
00044                         if ( p->str().find("ReferenceLine") != string::npos )
00045                         {
00046 
00047                                 cout <<"Can not define the Brezzi-Douglas-Marini element on a line"<<endl;
00048 
00049                         }
00050 
00051                         if ( p->str().find("ReferenceTetrahedron") != string::npos )
00052                         {
00053 
00054                                 cout <<"Can not define the Brezzi-Douglas-Marini element on a tetrahedron, see Nedelec2Hdiv"<<endl;
00055 
00056                         }
00057                         else if ( p->str().find("Triangle") != string::npos )
00058                         {
00059 
00060                                 description = istr("BrezziDouglasMarini_", order) + "_2D";
00061 
00062                                 Triangle& triangle = (Triangle&)(*p);
00063                                 GiNaC::lst equations;
00064                                 GiNaC::lst variables;
00065                                 GiNaC::lst polynom_space = bernsteinv(2,order, triangle, "b");
00066                                 GiNaC::ex pspace = polynom_space.op(0);
00067 
00068                                 variables = collapse(GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)));
00069 
00070                                 GiNaC::symbol t("t");
00071                                 GiNaC::ex dofi;
00072                                 // dofs related to edges
00073                                 for (int i=0; i< 3; i++)
00074                                 {
00075                                         Line line = triangle.line(i);
00076                                         GiNaC::lst normal_vec = normal(triangle, i);
00077                                         GiNaC::ex Vn = inner(pspace, normal_vec);
00078                                         GiNaC::lst points = interior_coordinates(line, order);
00079 
00080                                         GiNaC::ex point;
00081                                         for (unsigned int j=0; j< points.nops(); j++)
00082                                         {
00083                                                 point = points.op(j);
00084                                                 dofi = Vn.subs(x == point.op(0)).subs(y == point.op(1));
00085                                                 GiNaC::ex eq = dofi == GiNaC::numeric(0);
00086                                                 equations.append(eq);
00087                                                 dofs.insert(dofs.end(), GiNaC::lst(point,nsymb));
00088                                         }
00089                                 }
00090 
00091 //                              std::cout <<"no variables "<<variables.nops()<<std::endl;
00092 //                              std::cout <<"no equations "<<equations.nops()<<std::endl;
00093 
00094                                 // invert the matrix:
00095                                 // GiNaC has a bit strange way to invert a matrix.
00096                                 // It solves the system AA^{-1} = Id.
00097                                 // It seems that this way is the only way to do
00098                                 // properly with the solve_algo::gauss flag.
00099                                 //
00100                                 GiNaC::matrix b; GiNaC::matrix A;
00101                                 matrix_from_equations(equations, variables, A, b);
00102 
00103                                 unsigned int ncols = A.cols();
00104                                 GiNaC::matrix vars_sq(ncols, ncols);
00105 
00106                                 // matrix of symbols
00107                                 for (unsigned r=0; r<ncols; ++r)
00108                                         for (unsigned c=0; c<ncols; ++c)
00109                                                 vars_sq(r, c) = GiNaC::symbol();
00110 
00111                                 GiNaC::matrix id(ncols, ncols);
00112 
00113                                 // identity
00114                                 const GiNaC::ex _ex1(1);
00115                                 for (unsigned i=0; i<ncols; ++i)
00116                                         id(i, i) = _ex1;
00117 
00118                                 // invert the matrix
00119                                 GiNaC::matrix m_inv(ncols, ncols);
00120                                 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
00121 
00122                                 for (unsigned int i=0; i<dofs.size(); i++)
00123                                 {
00124                                         b.let_op(i) = GiNaC::numeric(1);
00125                                         GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
00126 
00127                                         GiNaC::lst subs;
00128                                         for (unsigned int ii=0; ii<xx.nops(); ii++)
00129                                         {
00130                                                 subs.append(variables.op(ii) == xx.op(ii));
00131                                         }
00132                                         GiNaC::ex Nj1 = pspace.op(0).subs(subs);
00133                                         GiNaC::ex Nj2 = pspace.op(1).subs(subs);
00134                                         Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
00135                                         b.let_op(i) = GiNaC::numeric(0);
00136                                 }
00137                         }
00138                 }
00139 
00140         }
00141 
00142 }                                                                // namespace SyFi
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines