/* $Id: FEMSolution.java,v 1.1 2000/03/14 14:39:09 thomasj Exp $ */ /****************************************************************************** * FEMSolution -- contains the finite elements solver. * @see FluidSolution * @see Solution ******************************************************************************/ public class FEMSolution extends FluidSolution{ /** Creates a FEMSolution object. @param mesh The mesh of the solution */ public FEMSolution(Mesh mesh){ super(mesh); } // FEMSolution //TAG_FEMInit_TAG// /** Discretize the initial Shape function and initialize the moments @param The initial shape to be approximated ****************************************************************************/ public void discretize(ShapeFunction function){ for (int j = 0; j < mesh_.size(); j++) //Discretizing is very simple f[j] = function.getValue(mesh_, j); // with normalized FEM initialMoments[0]=calculateMoments(0); //Set conserved quantities initialMoments[1]=calculateMoments(1); } // discretize //======================================================================= // Project 1:2 //======================================================================= /** Quadrature parameters */ int N_quad; double[] x_quad; double[] weigth; /** Initialize Gaussian Quadrature on the interval [0,1]*/ private void GaussianQuad() { N_quad=2; x_quad = new double[N_quad]; weigth = new double[N_quad]; x_quad[0]=0.5*(1-Math.sqrt(1./3.)); x_quad[1]=0.5*(1+Math.sqrt(1./3.)); weigth[0]=1; weigth[1]=1; } /** Evaluate the FE-basis function, number "i" at position X. */ private double FEbasis(int i,double X) { if (X-x[i]>0) { if (i==0 && X-x[f.length-1]>0) return ( X-x[f.length-1] )/dx[f.length-1]; return 1-(X-x[i])/dx[i]; } return ( X-x[i-1] )/dx[i-1]; } /** Evaluate the derivative of the FE-basis function, number "i" at position X. */ private double FEbasisPrim(int i,double X) { if (X-x[i]<0) return 1/dx[i-1]; if (i==0 && X-x[f.length-1]>0 ) return 1/dx[f.length-1]; return -1/dx[i]; } /** The diffusion coefficient at position X. */ private double Diffusion(double X, double diffusCo) { return diffusCo; } /** The velocity at position X. */ private double Velocity(double X, double velocity) { return velocity; } //====== End Project 1:2 =================================== /** Advance the solution forward one step in time. The equation is set by the internal variable equation_ and the numerical scheme by the internal variable scheme_ @param runData List of run parameters @see Data#TIMSTP @see Data#VELCTY @see Data#DIFFUS @see Data#THETA @see Data#TUNINT @return True if a scheme is implemented for that equation ****************************************************************************/ public boolean next(Data runData) { double timeStep = runData.getDouble(Data.TIMSTP); //Parameters double velocity = runData.getDouble(Data.VELCTY); double diffusCo = runData.getDouble(Data.DIFFUS); double theta = runData.getDouble(Data.THETA); double tune = runData.getDouble(Data.TUNINT); double alpha=timeStep*diffusCo/(dx[0]*dx[0]); //These are only constant if double beta =timeStep*velocity/(dx[0] ); // the problem and mesh are // homogeneous int n=f.length-1; boolean isDefined=false; if(equation_.equals(Data.ADVECTION) && scheme_.equals(Data.FEMMINE)){ //TAG_FEMAdvMine_TAG// // if(equation_.equals(Data.ADVECTION) // && scheme_.equals(Data.TUNFEM)){ //TAG_FEMAdv_TAG// //======================================================================= // FEM - Advection/diffusion, tunable integration t=1. FD // t=1/3 linear FEM // t=0. constant FEM //======================================================================= BandMatrix a = new BandMatrix(3, f.length); BandMatrix b = new BandMatrix(3, f.length); double[] c = new double[f.length]; double h = dx[0]; double htm = h*(1-tune)/4; double htp = h*(1+tune)/4; for (int i=0; i<=n; i++) { // Standard JBONE FEM-matrix. a.setL(i, htm +h*(-0.5*beta -alpha)* theta ); a.setD(i,2*(htp +h*( alpha)* theta )); a.setR(i, htm +h*( 0.5*beta -alpha)* theta ); b.setL(i, htm +h*(-0.5*beta -alpha)*(theta-1) ); b.setD(i,2*(htp +h*( alpha)*(theta-1))); b.setR(i, htm +h*( 0.5*beta -alpha)*(theta-1) ); } c=b.dot(f); //Right hand side c[0]=c[0]+b.getL(0)*f[n]; // with periodicity c[n]=c[n]+b.getR(n)*f[0]; fp=a.solve3(c); //Solve linear problem isDefined = false; } else if(equation_.equals(Data.ADVECTION) && scheme_.equals(Data.TUNFEM)){ //TAG_FEMAdv_TAG// // } else if(equation_.equals(Data.ADVECTION) // && scheme_.equals(Data.FEMMINE)){ //TAG_FEMAdvMine_TAG// //======================================================================= // FEM - Implement your own scheme here //======================================================================= BandMatrix a = new BandMatrix(3, f.length); BandMatrix b = new BandMatrix(3, f.length); double[] c = new double[f.length]; //======================================================================= // Project 2:2 //======================================================================= GaussianQuad(); int ip; double[] A_sub= new double[4]; double[] B_sub= new double[4]; double t0t=theta*timeStep; double t1t=(1-theta)*timeStep; a.setD(0,0); b.setD(0,0); for (int i=0; i<=n; i++) { // Initialize diagonal a.setD(i,0); b.setD(i,0); } for (int i=0; i<=n; i++) { // Loop over grid. if (i1) System.out.println(i+"i "+e_i+" "+xx); if (Math.abs(e_ip)>1) System.out.println(ip+"ip "+e_ip+" "+xx); double vel=Velocity(xx,velocity); double dif=Diffusion(xx,diffusCo); A_sub[0]+= e_i * e_i * weigth[j]*dx[i]/2; B_sub[0]+=(vel*e_i + dif*ep_i )*ep_i * weigth[j]*dx[i]/2; A_sub[1]+= e_i * e_ip * weigth[j]*dx[i]/2; B_sub[1]+=(vel*e_i + dif*ep_i )*ep_ip * weigth[j]*dx[i]/2; A_sub[2]+= e_ip * e_i * weigth[j]*dx[i]/2; B_sub[2]+=(vel*e_ip + dif*ep_ip)*ep_i * weigth[j]*dx[i]/2; A_sub[3]+= e_ip * e_ip * weigth[j]*dx[i]/2; B_sub[3]+=(vel*e_ip + dif*ep_ip)*ep_ip * weigth[j]*dx[i]/2; } a.setD(i , A_sub[0] + B_sub[0]*t0t + a.getD(i)); a.setR(i , A_sub[1] + B_sub[1]*t0t); a.setL(ip, A_sub[2] + B_sub[2]*t0t); a.setD(ip, A_sub[3] + B_sub[3]*t0t + a.getD(ip)); b.setD(i , A_sub[0] - B_sub[0]*t1t + b.getD(i)); b.setR(i , A_sub[1] - B_sub[1]*t1t); b.setL(ip, A_sub[2] - B_sub[2]*t1t); b.setD(ip, A_sub[3] - B_sub[3]*t1t + b.getD(ip)); } //====== End Project 2:2 =================================== c=b.dot(f); //Right hand side c[0]=c[0]+b.getL(0)*f[n]; // with periodicity c[n]=c[n]+b.getR(n)*f[0]; fp=a.solve3(c); //Solve linear problem isDefined = true; } // if if(isDefined) shiftLevels(); return isDefined; } // next /** Take the solution backward one step to initialize schemes with 3 time levels; not really appropriate in this context. @param runData List of run parameters @return False since it is not used here ****************************************************************************/ public boolean previous(Data runData){ return false; } } // class FEMSolution