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     // 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 

Generated on Mon Aug 31 16:16:45 2009 for SyFi by  doxygen 1.5.9