SyFi 0.3
|
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 "Bubble.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 Bubble::Bubble(): StandardFE() 00016 { 00017 description = "Bubble"; 00018 } 00019 00020 Bubble::Bubble(Polygon& p, unsigned int order): StandardFE(p, order) 00021 { 00022 compute_basis_functions(); 00023 } 00024 00025 void Bubble:: 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 ( p->str().find("ReferenceLine") != string::npos ) 00038 { 00039 Ns.insert(Ns.end(), x*(1-x)); 00040 description = "Bubble_1D"; 00041 } 00042 else if ( p->str().find("Triangle") != string::npos ) 00043 { 00044 00045 GiNaC::ex b = barycenter_triangle(p->vertex(0), p->vertex(1), p->vertex(2)); 00046 GiNaC::ex N = GiNaC::numeric(1); 00047 for (unsigned int d=0; d< b.nops(); d++) 00048 { 00049 N = N*b.op(d).rhs(); 00050 } 00051 Ns.insert(Ns.end(), N); 00052 00053 description = "Bubble_2D"; 00054 00055 } 00056 else if ( p->str().find("Tetrahedron") != string::npos ) 00057 { 00058 GiNaC::ex b = barycenter_tetrahedron(p->vertex(0), p->vertex(1), 00059 p->vertex(2), p->vertex(3)); 00060 GiNaC::ex N = GiNaC::numeric(1); 00061 for (unsigned int d=0; d< b.nops(); d++) 00062 { 00063 N = N*b.op(d).rhs(); 00064 } 00065 Ns.insert(Ns.end(), N); 00066 00067 description = "Bubble_3D"; 00068 } 00069 00070 // create and insert dof 00071 GiNaC::lst midpoint = GiNaC::lst(); 00072 for (unsigned int d=0; d< p->vertex(1).nops(); d++) 00073 { 00074 midpoint.append(GiNaC::numeric(0)); 00075 } 00076 for (unsigned int i=1; i<= p->no_vertices(); i++) 00077 { 00078 for (unsigned int d=0; d< p->vertex(i-1).nops(); d++) 00079 { 00080 midpoint.let_op(d) += p->vertex(i-1).op(d); 00081 } 00082 } 00083 00084 for (unsigned int d=0; d< p->vertex(1).nops(); d++) 00085 { 00086 midpoint.let_op(d) = midpoint.op(d)/p->no_vertices(); 00087 } 00088 00089 dofs.insert(dofs.end(), midpoint); 00090 } 00091 00092 } // namespace SyFi