BrezziDouglasMarini.cpp
Go to the documentation of this file.00001
00002
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
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
00092
00093
00094
00095
00096
00097
00098
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
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
00114 const GiNaC::ex _ex1(1);
00115 for (unsigned i=0; i<ncols; ++i)
00116 id(i, i) = _ex1;
00117
00118
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 }