arnoldfalkwinther.cpp
Go to the documentation of this file.00001 #include <SyFi.h>
00002 #include <fstream>
00003
00004
00005
00006 using namespace GiNaC;
00007 using namespace SyFi;
00008 using namespace std;
00009
00010
00011
00012 int main() {
00013
00014
00015 archive ar;
00016
00017 initSyFi(3);
00018 ReferenceTetrahedron tetrahedron;
00019 int order = 2;
00020
00021 ArnoldFalkWintherWeakSymSigma sigma_fe;
00022 sigma_fe.set_order(order);
00023 sigma_fe.set_polygon(tetrahedron);
00024 sigma_fe.compute_basis_functions();
00025 for (int i=0; i<sigma_fe.nbf(); i++) {
00026 cout <<"sigma_fe.N("<<i<<")="<<sigma_fe.N(i)<<endl; ;
00027 cout <<"div(sigma_fe.N("<<i<<"))="<<div(sigma_fe.N(i))<<endl; ;
00028 cout <<"sigma_fe.dof("<<i<<"))="<<sigma_fe.dof(i)<<endl; ;
00029
00030 ar.archive_ex(sigma_fe.N(i), istr("sN", i).c_str());
00031
00032 }
00033
00034 map<pair<int,int>,ex> A;
00035 pair<int,int> index;
00036 for (int i=0; i<sigma_fe.nbf(); i++) {
00037 index.first = i;
00038 for (int j=0; j<sigma_fe.nbf(); j++) {
00039 index.second = j;
00040 ex integrand = inner(sigma_fe.N(i), sigma_fe.N(j));
00041 A[index] = tetrahedron.integrate(integrand);
00042 cout <<"A["<<i<<","<<j<<"]="<<A[index]<<endl;
00043
00044 ar.archive_ex(A[index], istr("A", i,j).c_str());
00045 }
00046 }
00047
00048
00049
00050
00051 ArnoldFalkWintherWeakSymU u_fe;
00052 u_fe.set_order(order);
00053 u_fe.set_polygon(tetrahedron);
00054 u_fe.compute_basis_functions();
00055 for (int i=0; i<u_fe.nbf(); i++) {
00056 cout <<"u_fe.N("<<i<<")="<<u_fe.N(i)<<endl; ;
00057 ar.archive_ex(u_fe.N(i), istr("uN", i).c_str());
00058 }
00059
00060 map<pair<int,int>,ex> B;
00061 for (int i=0; i<sigma_fe.nbf(); i++) {
00062 index.first = i;
00063 for (int j=0; j<u_fe.nbf(); j++) {
00064 index.second = j;
00065 ex integrand = inner(div(sigma_fe.N(i)), u_fe.N(j));
00066 B[index] = tetrahedron.integrate(integrand);
00067 cout <<"B["<<i<<","<<j<<"]="<<B[index]<<endl;
00068 }
00069 }
00070
00071
00072
00073
00074 ArnoldFalkWintherWeakSymP p_fe;
00075 p_fe.set_order(order);
00076 p_fe.set_polygon(tetrahedron);
00077 p_fe.compute_basis_functions();
00078 for (int i=0; i<p_fe.nbf(); i++) {
00079 cout <<"p_fe.N("<<i<<")="<<p_fe.N(i)<<endl; ;
00080 ar.archive_ex(p_fe.N(i), istr("pN", i).c_str());
00081 }
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 return 0;
00094 }
00095