00001 #include <SyFi.h>
00002 #include <fstream>
00003
00004 using namespace GiNaC;
00005 using namespace SyFi;
00006 using namespace std;
00007
00008 void compute_poisson_element_matrix(
00009 FE& fe,
00010 Dof& dof,
00011 std::map<std::pair<unsigned int,unsigned int>, ex>& A)
00012 {
00013 std::pair<unsigned int,unsigned int> index;
00014 Polygon& domain = fe.get_polygon();
00015
00016 ex ujs = symbolic_matrix(1,fe.nbf(), "u");
00017 ex u;
00018 for (unsigned int k=0; k< fe.nbf(); k++) {
00019 u += ujs.op(k)*fe.N(k);
00020 }
00021
00022
00023 for (unsigned int i=0; i< fe.nbf() ; i++) {
00024 index.first = i;
00025 ex Fi = inner(grad(u), grad(fe.N(i)));
00026 for (unsigned int j=0; j< fe.nbf() ; j++) {
00027 index.second = j;
00028 symbol uj = ex_to<symbol>(ujs.op(j));
00029 ex nabla = Fi.diff(uj,1);
00030 ex Aij = domain.integrate(nabla);
00031 A[index] += Aij;
00032 }
00033 }
00034 }
00035
00036 void compute_nlconvdiff_element_matrix(
00037 FE& fe,
00038 Dof& dof,
00039 std::map<std::pair<unsigned int,unsigned int>, ex>& A)
00040 {
00041 std::pair<unsigned int,unsigned int> index;
00042 Polygon& domain = fe.get_polygon();
00043
00044
00045 for (unsigned int i=0; i< fe.nbf() ; i++) {
00046 dof.insert_dof(1,i,fe.dof(i));
00047 }
00048
00049
00050 ex UU = matrix(2,1,lst(0,0));
00051 ex ujs = symbolic_matrix(1,fe.nbf(), "u");
00052 for (unsigned int k=0; k< fe.nbf(); k++) {
00053 UU +=ujs.op(k)*fe.N(k);
00054 }
00055
00056
00057 matrix U = ex_to<matrix>(UU.evalm());
00058
00059 for (unsigned int i=0; i< fe.nbf() ; i++) {
00060 index.first = dof.glob_dof(fe.dof(i));
00061
00062
00063 ex gradU = grad(U);
00064 ex Fi_diffusion = inner(gradU, grad(fe.N(i)));
00065
00066
00067 ex Ut = U.transpose();
00068 ex UgradU = (Ut*gradU).evalm();
00069 ex Fi_convection = inner(UgradU, fe.N(i), true);
00070
00071
00072 ex Fi = Fi_convection + Fi_diffusion;
00073
00074
00075
00076
00077 for (unsigned int j=0; j< fe.nbf() ; j++) {
00078 index.second = dof.glob_dof(fe.dof(j));
00079 symbol uj = ex_to<symbol>(ujs.op(j));
00080 ex Jij = Fi.diff(uj,1);
00081 ex Aij = domain.integrate(Jij);
00082 A[index] += Aij;
00083 }
00084 }
00085 }
00086
00087
00088
00089
00090 int main() {
00091
00092 initSyFi(2);
00093
00094 Triangle T(lst(0,0), lst(1,0), lst(0,1), "t");
00095 int order = 2;
00096
00097
00098
00099
00100 Lagrange fe;
00101 fe.set_order(order);
00102 fe.set_polygon(T);
00103 fe.compute_basis_functions();
00104
00105 Dof dof1;
00106 std::map<std::pair<unsigned int,unsigned int>, ex> A1;
00107 compute_poisson_element_matrix(fe,dof1,A1);
00108 print(A1);
00109
00110
00111
00112 VectorLagrange vfe;
00113 vfe.set_order(order);
00114 vfe.set_size(2);
00115 vfe.set_polygon(T);
00116 vfe.compute_basis_functions();
00117 usage(vfe);
00118
00119 Dof dof2;
00120 std::map<std::pair<unsigned int,unsigned int>, ex> A2;
00121 compute_nlconvdiff_element_matrix(vfe,dof2,A2);
00122 cout <<"standard output"<<endl;
00123 print(A2);
00124 cout <<"LaTeX output"<<endl;
00125 cout <<latex;
00126 print(A2);
00127 cout <<"Python output"<<endl;
00128 cout <<python;
00129 print(A2);
00130 cout <<"C output"<<endl;
00131 cout <<csrc;
00132 print(A2);
00133
00134
00135
00136
00137
00138 archive ar;
00139 map<std::pair<unsigned int,unsigned int>,ex>::iterator iter;
00140 for (iter = A1.begin(); iter != A1.end() ; iter++) {
00141 ar.archive_ex((*iter).second, istr("A1_",
00142 (*iter).first.first,
00143 (*iter).first.second).c_str());
00144 }
00145 for (iter = A2.begin(); iter != A2.end() ; iter++) {
00146 ar.archive_ex((*iter).second, istr("A2_",
00147 (*iter).first.first,
00148 (*iter).first.second).c_str());
00149 }
00150
00151 ofstream vfile("nljacobian_ex.gar.v");
00152 vfile << ar; vfile.close();
00153 if(!compare_archives("nljacobian_ex.gar.v", "nljacobian_ex.gar.r")) {
00154 cerr << "Failure!" << endl;
00155 return -1;
00156 }
00157
00158 return 0;
00159 }
00160