/* $Id: FEMSolution.java,v 1.1 1999/12/09 14:59:34 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 */ private double V[]; //Potential for Schroedinger eq. 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){ if (scheme_.equals(Data.FEMMINE2)) { setPotential(); //Potential for Schrödinger eq. } 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 /** 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.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++) { 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 = true; } else if(equation_.equals(Data.ADVECTION) && scheme_.equals(Data.FEMMINE)){ //TAG_FEMAdvMine_TAG// //======================================================================= // FEM - Diffusion in a cylinder (x=r) //======================================================================= 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 g = timeStep*diffusCo; for (int i=0; i<=n; i++) { a.setL(i,h*h*1./12+h*1./6*i-(g*0.5+g/h*i)*theta); a.setD(i,h*2./3*i+g/h*2.*i*theta); a.setR(i,-h*h*1./12+h*i*1./6+(g*0.5-g/h*i)*theta); b.setL(i,h*h*1./12+h*1./6*i-(g*0.5+g/h*i)*(theta-1)); b.setD(i,h*2./3*i+g/h*2.*i*(theta-1)); b.setR(i,-h*h*1./12+h*i*1./6+(g*0.5-g/h*i)*(theta-1)); } c=b.dot(f); //Right hand side a.setL(0,0.); //Neumann boundary condition dT/dr=0, T0=T1, r=0 a.setD(0,1.); a.setR(0,-1); a.setL(n,-1.); //Neumann boundary condition dT/dr=0,Tn=Tn-1, r=rc a.setD(n,1.); a.setR(n,0.); c[0]=0; // with periodicity c[n]=0; fp=a.solve3(c); //Solve linear problem isDefined = true; } if(equation_.equals(Data.ADVECTION) && scheme_.equals(Data.FEMMINE2)){ //TAG_FEMAdv_TAG// //======================================================================= // FEM - Advection/diffusion, tunable integration t=1. FD // Schrödinger equation->Complex-values t=1/3 linear FEM // t=0. constant FEM //======================================================================= BandMatrixC a = new BandMatrixC(3, h.length); BandMatrixC b = new BandMatrixC(3, h.length); Complex[] c = new Complex[h.length]; Complex I = new Complex(0.,1.); double g = dx[0]; //Matrix elements double htm = g*(1-tune)/4; Complex chtm = new Complex(htm,0.); double htp = g*(1+tune)/4; double htp2 =2.* htp; Complex chtp2 = new Complex(htp2,0.); Complex ctmp = new Complex(); Complex ctmp1 = new Complex(); for (int i=0; i<=n; i++) { double tmpm = (timeStep*V[i]*htm - timeStep/g)*theta; Complex ctmpm = new Complex(tmpm,0.); double tmpp = (timeStep*V[i]*htp2 + 2*timeStep/g)* theta; Complex ctmpp = new Complex(tmpp,0.); double tmpm1 = (timeStep*V[i]*htm - timeStep/g)*(theta-1); Complex ctmpm1 = new Complex(tmpm1,0.); double tmpp1 = (timeStep*V[i]*htp2 + 2*timeStep/g)*(theta-1); Complex ctmpp1 = new Complex(tmpp1,0.); ctmp = chtm.mul(I); a.setL(i,ctmp.sub(ctmpm)); ctmp1 = chtp2.mul(I); a.setD(i,ctmp1.sub(ctmpp)); ctmp = chtm.mul(I); a.setR(i,ctmp.sub(ctmpm)); ctmp = chtm.mul(I); b.setL(i,ctmp.sub(ctmpm1)); ctmp1 = chtp2.mul(I); b.setD(i,ctmp1.sub(ctmpp1)); ctmp = chtm.mul(I); b.setR(i,ctmp.sub(ctmpm1)); } c=b.dot(h); //Right hand side Complex tmp = new Complex(); //with periodicity tmp=b.getL(0).mul(h[n]); c[0]=c[0].add(tmp); Complex tmp1 = new Complex(); tmp1=b.getR(n).mul(h[0]); c[n]=c[n].add(tmp1); hp=a.solve3(c); //Solve linear problem for(int j=0; j<=n; j++){ //To put the real value of h in f to plot fp[j]=hp[j].re(); //System.out.println("j="+j+" hp["+j+"]=" +hp[j]); //fp[j]=hp[j].abs()*hp[j].abs(); //The probability } // for(int j=0; j<=n; j++) fp[j]=hp[j].re(); // Test loop 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){ if (scheme_.equals(Data.FEMMINE2)) { for (int j = 0; j < mesh_.size(); j++) { double position = runData.getDouble(Data.IC1POS); //Wavepacket double width = runData.getDouble(Data.IC1WID); double wavelength = runData.getDouble(Data.IC1LEN); double amplitude = runData.getDouble(Data.IC1AMP); double norm=Math.sqrt(Math.PI/width); double xx0 = (mesh_.point(j) - position); Complex arg1 =new Complex(0.,2.*Math.PI/wavelength*xx0); double arg2 =-xx0*xx0/(4.*width); Complex psi = new Complex(arg1.exp().scale(norm*Math.exp(arg2)) ); h[j] = psi; //Change amplitude to norm? } //Real is plotted initialMoments[0]=calculateMoments(0); //Set conserved quantities initialMoments[1]=calculateMoments(1); } return true; } /** Distcretize the potential **/ public void setPotential(){ int n=f.length-1; V = new double[f.length]; //Create the object - potential array //----Barrier double steep =10.; // Range for tanh(x), x in [-steep;+steep] double T=2* Math.PI; // Constant double A=1.5*T; //Amplitude of the potential //double A=0.; // Amplitude eual to zero! //double A=3.0*T; //Large potential for(int k=0; k<=n; k++){ double dx = 2.*steep/((double)n); double x = (double)k*dx-steep; double tanh = (Math.exp(x)-Math.exp(-x))/(Math.exp(x)+Math.exp(-x)); V[k] = A*Math.sqrt(0.5*(1.+tanh)); } //---Harmonic-periodic /* double dx = 0.; //Common variables double x = 0.; double ampl=81.; double half=32.; double tmp=0.; for(int k=0; k<=n; k++){ dx=2.*half/((double)n); x=0.5*((double)k*dx-half)/half; V[k]=ampl*(x*x); } */ } // class FEMSolution }