ElementComputations.cpp
Go to the documentation of this file.00001
00002
00003
00004 #include "ElementComputations.h"
00005 #include "tools.h"
00006
00007 using std::cout;
00008 using std::endl;
00009
00010 namespace SyFi
00011 {
00012
00013 void usage(FE& fe)
00014 {
00015 for (unsigned int i=0; i< fe.nbf(); i++)
00016 {
00017 cout <<"fe.N("<<i<<") = "<<fe.N(i)<<endl;
00018 cout <<"grad(fe.N("<<i<<")) = "<<grad(fe.N(i))<<endl;
00019 cout <<"fe.dof("<<i<<") = "<<fe.dof(i)<<endl;
00020 }
00021 }
00022
00023 void usage(FE& v_fe, FE& p_fe)
00024 {
00025 for (unsigned int i=0; i< v_fe.nbf(); i++)
00026 {
00027 cout <<"v_fe.N("<<i<<") = "<<v_fe.N(i)<<endl;
00028 cout <<"grad(v_fe.N("<<i<<")) = "<<grad(v_fe.N(i))<<endl;
00029 cout <<"v_fe.dof("<<i<<") = "<<v_fe.dof(i)<<endl;
00030 }
00031 for (unsigned int i=0; i< p_fe.nbf(); i++)
00032 {
00033 cout <<"p_fe.N("<<i<<")= "<<p_fe.N(i)<<endl;
00034 cout <<"p_fe.dof("<<i<<")= "<<p_fe.dof(i)<<endl;
00035 }
00036 }
00037
00038 void compute_Poisson_element_matrix(
00039 FE& fe,
00040 Dof& dof,
00041 std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A)
00042 {
00043 std::pair<unsigned int,unsigned int> index;
00044
00045
00046 for (unsigned int i=0; i< fe.nbf(); i++)
00047 {
00048 dof.insert_dof(1,i,fe.dof(i));
00049 }
00050
00051 Polygon& domain = fe.get_polygon();
00052
00053
00054 for (unsigned int i=0; i< fe.nbf(); i++)
00055 {
00056
00057 index.first = dof.glob_dof(fe.dof(i));
00058 for (unsigned int j=0; j< fe.nbf(); j++)
00059 {
00060
00061 index.second = dof.glob_dof(fe.dof(j));
00062
00063 GiNaC::ex nabla = inner(grad(fe.N(i)),
00064 grad(fe.N(j)));
00065
00066 GiNaC::ex Aij = domain.integrate(nabla);
00067 A[index] += Aij;
00068 }
00069 }
00070 }
00071
00072 void compute_Stokes_element_matrix(
00073 FE& v_fe,
00074 FE& p_fe,
00075 Dof& dof,
00076 std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A)
00077 {
00078 std::pair<unsigned int,unsigned int> index;
00079 std::pair<unsigned int,unsigned int> index2;
00080
00081
00082
00083 Polygon& domain = v_fe.get_polygon();
00084
00085
00086 for (unsigned int i=0; i< v_fe.nbf(); i++)
00087 {
00088 dof.insert_dof(1,i,v_fe.dof(i));
00089 }
00090 for (unsigned int i=0; i< p_fe.nbf(); i++)
00091 {
00092 dof.insert_dof(1,v_fe.nbf()+i,p_fe.dof(i));
00093 }
00094
00095
00096 for (unsigned int i=0; i< v_fe.nbf(); i++)
00097 {
00098
00099 index.first = dof.glob_dof(v_fe.dof(i));
00100 for (unsigned int j=0; j< v_fe.nbf(); j++)
00101 {
00102
00103 index.second = dof.glob_dof(v_fe.dof(j));
00104 GiNaC::ex nabla = inner(grad(v_fe.N(i)),
00105 grad(v_fe.N(j)));
00106
00107 GiNaC::ex Aij = domain.integrate(nabla);
00108 A[index] += Aij;
00109 }
00110 }
00111
00112
00113 for (unsigned int i=0; i< p_fe.nbf(); i++)
00114 {
00115
00116 index.first = dof.glob_dof(p_fe.dof(i));
00117 for (unsigned int j=0; j< v_fe.nbf(); j++)
00118 {
00119
00120 index.second=dof.glob_dof(v_fe.dof(j));
00121
00122 GiNaC::ex divV= -p_fe.N(i)*div(v_fe.N(j));
00123
00124 GiNaC::ex Aij = domain.integrate(divV);
00125 A[index] += Aij;
00126
00127
00128
00129 index2.first = index.second;
00130 index2.second = index.first;
00131 A[index2] += Aij;
00132 }
00133 }
00134 }
00135
00136 void compute_mixed_Poisson_element_matrix(
00137 FE& v_fe,
00138 FE& p_fe,
00139 Dof& dof,
00140 std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A)
00141 {
00142 std::pair<unsigned int,unsigned int> index;
00143 std::pair<unsigned int,unsigned int> index2;
00144
00145
00146
00147 Polygon& domain = v_fe.get_polygon();
00148
00149
00150 for (unsigned int i=0; i< v_fe.nbf(); i++)
00151 {
00152 dof.insert_dof(1,i,v_fe.dof(i));
00153 }
00154 for (unsigned int i=0; i< p_fe.nbf(); i++)
00155 {
00156 dof.insert_dof(1,v_fe.nbf()+i+1,p_fe.dof(i));
00157 }
00158
00159
00160 for (unsigned int i=0; i< v_fe.nbf(); i++)
00161 {
00162
00163 index.first = dof.glob_dof(v_fe.dof(i));
00164 for (unsigned int j=0; j< v_fe.nbf(); j++)
00165 {
00166
00167 index.second = dof.glob_dof(v_fe.dof(j));
00168
00169 GiNaC::ex mass = inner(v_fe.N(i),v_fe.N(j));
00170
00171 GiNaC::ex Aij = domain.integrate(mass);
00172 A[index] += Aij;
00173 }
00174 }
00175
00176
00177 for (unsigned int i=0; i< p_fe.nbf(); i++)
00178 {
00179
00180 index.first = dof.glob_dof(p_fe.dof(i));
00181 for (unsigned int j=0; j< v_fe.nbf(); j++)
00182 {
00183
00184 index.second=dof.glob_dof(v_fe.dof(j));
00185
00186 GiNaC::ex divV= -p_fe.N(i)*div(v_fe.N(j));
00187
00188 GiNaC::ex Aij = domain.integrate(divV);
00189 A[index] += Aij;
00190
00191
00192
00193 index2.first = index.second;
00194 index2.second = index.first;
00195 A[index2] += Aij;
00196 }
00197 }
00198 }
00199
00200 }