SyFi 0.3
|
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 // the polynomial spaces on the form: 00020 // first item: a0 + a1*x + a2*y + a3*x^2 + a4*x*y ... the polynom 00021 // second item: a0, a1, a2, ... the coefficents 00022 // third item 1, x, y, x^2 the basis 00023 // Could also do: 00024 // GiNaC::ex polynom_space = bernstein(order, t, "a"); 00025 00026 polynom = polynom_space.op(0); 00027 variables = ex_to<lst>(polynom_space.op(1)); 00028 // the variables a0,a1,a2 .. 00029 00030 ex Nj; 00031 // The bezier ordinates (in which the basis function should be either 0 or 1) 00032 lst points = bezier_ordinates(t,order); 00033 00034 // Loop over all basis functions Nj and all points. 00035 // Each basis function Nj is determined by a set of linear equations: 00036 // Nj(xi) = dirac(i,j) 00037 // This system of equations is then solved by lsolve 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 // The point xi 00043 ex point = points.op(i-1); 00044 // The equation Nj(x) = dirac(i,j) 00045 ex eq = polynom == dirac(i,j); 00046 // Substitute x = xi and y = yi and appended the equation 00047 // to the list of equations 00048 equations.append(eq.subs(lst(x == point.op(0) , y == point.op(1)))); 00049 } 00050 00051 00052 // We solve the linear system 00053 // cout <<"equations "<<equations<<endl; 00054 // cout <<"variables "<<variables<<endl; 00055 ex subs = lsolve(equations, variables); 00056 // Substitute to get the N_j 00057 Nj = polynom.subs(subs); 00058 basis_functions.append(Nj); 00059 cout <<"Nj "<<Nj<<endl; 00060 } 00061 00062 00063 00064 // regression test 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