yuemingl / Futureye_JIT

package edu.uta.futureye.tutorial;

import static edu.uta.futureye.function.FMath.C0;
import static edu.uta.futureye.function.FMath.grad;
import static edu.uta.futureye.function.FMath.x;
import static edu.uta.futureye.function.FMath.y;

import java.util.HashMap;

import edu.uta.futureye.algebra.intf.Matrix;
import edu.uta.futureye.algebra.intf.Vector;
import edu.uta.futureye.algebra.solver.external.SolverJBLAS;
import edu.uta.futureye.core.Mesh;
import edu.uta.futureye.core.NodeType;
import edu.uta.futureye.function.intf.MathFunc;
import edu.uta.futureye.io.MeshReader;
import edu.uta.futureye.io.MeshWriter;
import edu.uta.futureye.lib.assembler.Assembler;
import edu.uta.futureye.lib.element.FELinearTriangle;
import edu.uta.futureye.lib.weakform.WeakForm;
import edu.uta.futureye.util.Utils;

public class LaplaceGlobalAssemble {
	public void run() {
		// 1.Read in mesh
		MeshReader reader = new MeshReader("grids/triangle.grd");
		Mesh mesh = reader.read2DMesh();
		// Compute geometry relationship between nodes and elements
		mesh.computeNodeBelongsToElements();

		// 2.Mark boundary type(s)
		HashMap<NodeType, MathFunc> mapNTF = new HashMap<NodeType, MathFunc>();
		mapNTF.put(NodeType.Dirichlet, null);  //null => mark all boundary nodes
		mesh.markBorderNode(mapNTF);

		// 3. Weak form definition
		FiniteElement fe = new FELinearTriangle(); // Linear triangular finite element
		final MathFunc f = -2 * (x * x + y * y) + 36; //Right hand side (RHS)
		WeakForm wf = new WeakForm(fe,
				(u,v) -> grad(u, "x", "y").dot(grad(v, "x", "y")), 
				v -> f * v
			);
		wf.compile();

		// 4. Assembly and boundary condition(s)
		Assembler assembler = new Assembler(mesh, wf);
		assembler.assembleGlobal();
		Matrix stiff = assembler.getGlobalStiffMatrix();
		Vector load = assembler.getGlobalLoadVector();
		// Apply zero Dirichlet boundary condition
		Utils.imposeDirichletCondition(stiff, load, fe, mesh, C0);

		// 5. Solve the linear system
		SolverJBLAS solver = new SolverJBLAS();
		Vector u = solver.solveDGESV(stiff, load);
		System.out.println("u=");
		for (int i = 1; i <= u.getDim(); i++)
			System.out.println(String.format("%.3f ", u.get(i)));

		// 6. Output the result to a file with Techplot format
		MeshWriter writer = new MeshWriter(mesh);
		writer.writeTechplot("./tutorial/Laplace2D.dat", u);
	}

	public static void main(String[] args) {
		LaplaceGlobalAssemble ex1 = new LaplaceGlobalAssemble();
		ex1.run();
	}
}

MathFunc upwindCoef = new UserDefFunc() {
	//@Override
	public double apply(AssembleParam ap, double... args) {
		Node node1 = ap.element.nodes.at(ap.testDOFIdx);
		Node node2 = ap.element.nodes.at(ap.testDOFIdx==1?2:1);
		double coord1 = node1.coord(1);
		double coord2 = node2.coord(1);
		double upwindWeight = 0.0;
		if((coord2-coord1)*u > 0) {
			upwindWeight = -0.1;
		} else {
			upwindWeight = 0.1;
		}
		return upwindWeight;
	}
};
// Weak form definition
WeakForm wf = new WeakForm(
		new FELinearLine1D(),
		(c, v) -> k * c.diff("x") * v.diff("x") + u * c.diff("x") * (v + upwindCoef),
		(v)    -> C0
		);
wf.compile();

/**
 * Stokes problem is defined as
 *  -\nabla\cdot(k\nabla{\vec{u}}) + \nabla{p} = \vec{f}\\
 *  div~{\vec{u}} = 0
 *
 *  or written explicitly in each dimension:
 *   -k*(u1_xx+u1_yy) + p_x = f1
 *   -k*(u2_xx+u2_yy) + p_y = f2
 *   u1_x+u2_y              = 0
 * The weak form is
 *   find \vec{u} \in H_0^1(div;\Omega), p \in L_2(\Omega)
 *   such that, for all \vec{v} \in H_0^1(div;\Omega), q \in L_2(\Omega)
 *   
 *   (\nabla{\vec{v}},k\nabla{\vec{u}}) - (div~{\vec{v}},p) 
 *                   + (q,div~{\vec{u}}) = (\vec{v},\vec{f})
 *   or written explicitly:
 *   (v1_x,k*u1_x) + (v1_y,k*u1_y) + (v2_x,k*u2_x) + (v2_y,k*u2_y) 
 *                 - (v1_x+v2_y,p) + (q,u1_x+u2_y) = (v1*f1+v2*f2)
 * where
 *   \vec{u}=(u1,u2): velocity vector field
 *   \vec{f}=(f1,f2): body force
 *
 */
// Weak form definition
FEQuadraticV_LinearP fe = new FEQuadraticV_LinearP();
MathFunc k = C1;
VecMathFunc f = new SpaceVectorFunction(C0, C0);
VecWeakForm wf = new VecWeakForm(fe,
		(u, v) -> k * grad(u[1],"x","y" ).dot(grad(v[1],"x","y")) //   (v1_x,k*u1_x) + (v1_y,k*u1_y)
				+ k * grad(u[2],"x","y" ).dot(grad(v[2],"x","y")) // + (v2_x,k*u2_x) + (v2_y,k*u2_y) 
				- (v[1].diff("x")+v[2].diff("y"))*u[3]            // - (v1_x+v2_y,p) //where p=u[3]
				+ v[3]*(u[1].diff("x")+u[2].diff("y")),           // + (q,u1_x+u2_y) //where q=v[3]
		(v)-> v[1]*f[1] + v[2]*f[2]);
wf.compile();