fe_ex3.cpp
Go to the documentation of this file.00001 #include <SyFi.h>
00002 #include <fstream>
00003
00004 using namespace GiNaC;
00005 using namespace SyFi;
00006 using namespace std;
00007
00008 int main() {
00009
00010 initSyFi(2);
00011
00012 Triangle t(lst(0,0), lst(1,0), lst(0,1));
00013 int order = 2;
00014 ex polynom;
00015 lst variables;
00016 lst basis_functions;
00017
00018 ex polynom_space = pol(order, 2, "a");
00019
00020
00021
00022
00023
00024
00025
00026 polynom = polynom_space.op(0);
00027 variables = ex_to<lst>(polynom_space.op(1));
00028
00029
00030 ex Nj;
00031
00032 lst points = bezier_ordinates(t,order);
00033
00034
00035
00036
00037
00038 for (int j=1; j <= points.nops(); j++) {
00039 lst equations;
00040 int i=0;
00041 for (int i=1; i<= points.nops() ; i++ ) {
00042
00043 ex point = points.op(i-1);
00044
00045 ex eq = polynom == dirac(i,j);
00046
00047
00048 equations.append(eq.subs(lst(x == point.op(0) , y == point.op(1))));
00049 }
00050
00051
00052
00053
00054
00055 ex subs = lsolve(equations, variables);
00056
00057 Nj = polynom.subs(subs);
00058 basis_functions.append(Nj);
00059 cout <<"Nj "<<Nj<<endl;
00060 }
00061
00062
00063
00064
00065
00066 archive ar;
00067 for (int i=0; i< basis_functions.nops(); i++) {
00068 ar.archive_ex(basis_functions[i], istr("N",i).c_str());
00069 }
00070
00071 ofstream vfile("fe_ex3.gar.v");
00072 vfile << ar; vfile.close();
00073 if(!compare_archives("fe_ex3.gar.v", "fe_ex3.gar.r")) {
00074 cerr << "Failure!" << endl;
00075 return -1;
00076 }
00077
00078 return 0;
00079
00080 }
00081
00082