gradient calculated.

it outputs all gradients. its speed is fine without OPENMP.	

Signed-off-by: Yimin <[email protected]>

src/Solver/Solver.m

@@ -36,6 +36,10 @@ function delete(this)
 	[DX, DY] = Solver_('reference', this.id_, points);
 end
 
+function [gradX, gradY] = grad(this, soln, nodes, elems, DX, DY)
+	[gradX, gradY] = Solver_('grad', this.id_, soln, nodes, elems, DX, DY);
+end
+
 function [x] = solve(this, A, b)
 
 if (this.ilu == 1)

src/Solver/private/Solver.cc

@@ -89,6 +89,67 @@ void Solver::Reference(MatlabPtr& DX, MatlabPtr& DY, MatlabPtr Points){
 
 }
 
+void Solver::Gradient(MatlabPtr& GradX, MatlabPtr& GradY, MatlabPtr Solution, MatlabPtr Nodes,
+		MatlabPtr Elems, MatlabPtr DX, MatlabPtr DY){
+
+	auto numberofelem        = mxGetN(Elems);
+	auto numberofnodeperelem = mxGetM(Elems);
+
+	auto gradX_ptr = mxGetPr(GradX);
+	auto gradY_ptr = mxGetPr(GradY);
+
+	auto elem_ptr = (int32_t*)mxGetPr(Elems);
+	auto node_ptr = mxGetPr(Nodes);
+
+	auto sln_ptr  = mxGetPr(Solution);
+
+	auto DX_ptr   = mxGetPr(DX);
+	auto DY_ptr   = mxGetPr(DY);
+
+
+	Real_t temp_ptr[numberofnodeperelem];
+	Real_t gx_ptr[numberofnodeperelem];
+	Real_t gy_ptr[numberofnodeperelem];
+
+	size_t vertex_1, vertex_2, vertex_3;
+	Real_t det;
+	Real_t Jacobian[2][2];
+
+	for (size_t i = 0; i < numberofelem; i++) {
+
+		for (size_t j = 0; j < numberofnodeperelem; j++) {
+			temp_ptr[j] = sln_ptr[elem_ptr[i * numberofnodeperelem + j] - 1];
+		}
+
+
+		for (size_t j = 0; j < numberofnodeperelem; j++) {
+			gx_ptr[j] = 0.;
+			gy_ptr[j] = 0.;
+			for (size_t k = 0; k < numberofnodeperelem; k++) {
+				gx_ptr[j] += temp_ptr[k] * DX_ptr[j * numberofnodeperelem + k];
+				gy_ptr[j] += temp_ptr[k] * DY_ptr[j * numberofnodeperelem + k];
+			}
+		}
+		vertex_1 = elem_ptr[numberofnodeperelem*i] - 1;
+		vertex_2 = elem_ptr[numberofnodeperelem*i + 1] - 1;
+		vertex_3 = elem_ptr[numberofnodeperelem*i + 2] - 1;
+
+		Jacobian[0][0] = node_ptr[2*vertex_3 + 1] - node_ptr[2*vertex_1 + 1];
+		Jacobian[1][1] = node_ptr[2*vertex_2    ] - node_ptr[2*vertex_1    ];
+		Jacobian[0][1] = node_ptr[2*vertex_1 + 1] - node_ptr[2*vertex_2 + 1];
+		Jacobian[1][0] = node_ptr[2*vertex_1    ] - node_ptr[2*vertex_3    ];
+
+		// Orientation corrected.
+		det = Jacobian[0][0] * Jacobian[1][1] - Jacobian[0][1] * Jacobian[1][0];
+
+		for (size_t j = 0; j < numberofnodeperelem; j++) {
+			gradX_ptr[i * numberofnodeperelem + j] =
+					(Jacobian[0][0] * gx_ptr[j] + Jacobian[0][1] * gy_ptr[j])/det;
+			gradY_ptr[i * numberofnodeperelem + j] =
+					(Jacobian[1][0] * gx_ptr[j] + Jacobian[1][1] * gy_ptr[j])/det;
+		}
+	}
+}
 } /* namespace MEX */
 
 using namespace MEX;
@@ -110,6 +171,24 @@ MEX_DEFINE(delete) (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
 	Session<Solver>::destroy(input.get(0));
 }
 
+
+MEX_DEFINE(grad) (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
+	InputArguments input(nrhs,  prhs, 6);
+	OutputArguments output(nlhs, plhs, 2);
+
+	Solver* solver = Session<Solver>::get(input.get(0));
+
+	size_t numberofelem        = mxGetN(prhs[3]);
+	size_t numberofnodeperelem = mxGetM(prhs[3]);
+
+	plhs[0] = mxCreateNumericMatrix(numberofnodeperelem, numberofelem, mxDOUBLE_CLASS, mxREAL);
+	plhs[1] = mxCreateNumericMatrix(numberofnodeperelem, numberofelem, mxDOUBLE_CLASS, mxREAL);
+
+	solver->Gradient(plhs[0], plhs[1], CAST(prhs[1]),
+			CAST(prhs[2]),CAST(prhs[3]),
+			CAST(prhs[4]),CAST(prhs[5]));
+}
+
 MEX_DEFINE(reference) (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
 	InputArguments input(nrhs,  prhs, 2);
 	OutputArguments output(nlhs, plhs, 2);

src/Solver/private/Solver.h

@@ -34,8 +34,14 @@ class Solver {
 	/*
 	 * public methods
 	 */
+	// return reference
 	void Reference(MatlabPtr &DX, MatlabPtr &DY,
 			MatlabPtr Points);
+	// return gradient on each element
+	void Gradient(MatlabPtr& GradX, MatlabPtr& GradY,MatlabPtr Solution, MatlabPtr Nodes, MatlabPtr Elems,
+			MatlabPtr DX, MatlabPtr DY);
+	// return Neumann data, unfinished yet
+	void Neumann();
 };
 
 } /* namespace MEX */