Improvement in shape function computation

FEM/CMakeLists.txt

@@ -49,6 +49,7 @@ set( HEADERS
 	tools.h
 	vtk.h
 	OutputTools.h
+	ShapeFunctionPool.h
 )
 
 set( SOURCES
@@ -102,6 +103,7 @@ set( SOURCES
 	Stiff_Bulirsch-Stoer.cpp
 	tools.cpp
 	vtk.cpp
+	ShapeFunctionPool.cpp
 )
 
 if(NOT OGS_LSOLVER STREQUAL PETSC)

FEM/DUMUX.cpp

@@ -1276,7 +1276,7 @@ void CDUMUXData::WriteDataToGeoSys(CRFProcess* m_pcs)
 			m_ele = m_msh->ele_vector[i]; // get element
 			if (m_ele->GetMark()) // Marked for use
 			{ // Configure Element for interpolation of node velocities to GP velocities
-				fem->ConfigElement(m_ele, m_pcs->m_num->ele_gauss_points);
+				fem->ConfigElement(m_ele);
 
 				std::string tempstring;
 				tempstring = "";

FEM/Output.cpp

@@ -18,6 +18,7 @@
 #include <string>
 
 #include <cfloat> // DBL_EPSILON
+#include <algorithm> // remove-if
 
 #include "BuildInfo.h"
 #include "FEMIO/GeoIO.h"
@@ -477,6 +478,7 @@ ios::pos_type COutput::Read(std::ifstream& in_str, const GEOLIB::GEOObjects& geo
 				}
 				if (Keyword(line_string))
 					return position;
+				std::remove_if(line_string.begin(), line_string.end(), ::isspace);
 				mfp_value_vector.push_back(line_string);
 			}
 
@@ -3722,13 +3724,12 @@ void COutput::CalculateTotalFlux(CFEMesh* msh, vector<long>& nodes_on_geo, vecto
 	double fac, nodesFVal[8], nodesFVal_adv[8], flux[3]; // , poro;
 	// CMediumProperties *MediaProp;
 
-	int Axisymm = 1; // ani-axisymmetry
-	if (msh->isAxisymmetry())
-		Axisymm = -1; // Axisymmetry is true
-
 	CElem* elem = NULL;
 	CElem* face = new CElem(1);
-	FiniteElement::CElement* element = new FiniteElement::CElement(Axisymm * msh->GetCoordinateFlag());
+
+	FiniteElement::CElement* fem_assembler = m_pcs_flow->getLinearFEMAssembler();
+	assert(fem_assembler);
+
 	CNode* e_node = NULL;
 	CElem* e_nei = NULL;
 	set<long> set_nodes_on_geo;
@@ -3796,11 +3797,12 @@ void COutput::CalculateTotalFlux(CFEMesh* msh, vector<long>& nodes_on_geo, vecto
 				fac = 0.5; // Not a surface face
 			face->SetFace(elem, j);
 			face->SetOrder(msh->getOrder());
+			face->FillTransformMatrix();
 			face->ComputeVolume();
-			face->SetNormalVector();
+			face->SetNormalVector();  // to get it directly from TransformMatrix
 			face->DirectNormalVector();
-			element->setOrder(msh->getOrder() + 1);
-			element->ConfigElement(face, m_pcs_flow->m_num->ele_gauss_points, true); // 2D fem
+			fem_assembler->setOrder(msh->getOrder() + 1);
+			fem_assembler->ConfigElement(face, true); // 2D fem
 
 			for (k = 0; k < nfn; k++)
 			{
@@ -3818,8 +3820,8 @@ void COutput::CalculateTotalFlux(CFEMesh* msh, vector<long>& nodes_on_geo, vecto
 					                   * mfp_vector[0]->SpecificHeatCapacity() * mfp_vector[0]->Density();
 			}
 			///
-			element->FaceNormalFluxIntegration(elements_at_geo[i], nodesFVal, nodesFVal_adv, nodesFace, face, m_pcs,
-			                                   face->normal_vector);
+			fem_assembler->FaceNormalFluxIntegration(elements_at_geo[i], nodesFVal, nodesFVal_adv, nodesFace, face,
+			                                         m_pcs, face->normal_vector);
 			for (k = 0; k < nfn; k++)
 			{
 				e_node = elem->GetNode(nodesFace[k]);
@@ -3841,7 +3843,6 @@ void COutput::CalculateTotalFlux(CFEMesh* msh, vector<long>& nodes_on_geo, vecto
 	NVal_diff.clear();
 	NVal_adv.clear();
 	G2L.clear();
-	delete element;
 	delete face;
 }
 

FEM/OutputTools.h

@@ -106,13 +106,11 @@ inline double getElementMMP(int mmp_id, MeshLib::CElem* ele, CRFProcess* m_pcs)
 {
 	double gp[3] = {.0, .0, .0};
 	double theta = 1.0;
-	int gp_r, gp_s, gp_t;
 	ele->SetOrder(false);
 	CFiniteElementStd* fem = m_pcs->GetAssember();
-	fem->ConfigElement(ele, m_pcs->m_num->ele_gauss_points, false);
+	fem->ConfigElement(ele, false);
 	fem->Config();
-	fem->SetGaussPoint(0, gp_r, gp_s, gp_t);
-	fem->ComputeShapefct(1);
+	fem->getShapeFunctionCentroid();
 	CMediumProperties* mmp = mmp_vector[ele->GetPatchIndex()];
 	double val = ELEMENT_MMP_VALUES::getValue(mmp, mmp_id, ele->GetIndex(), gp, theta);
 	return val;

FEM/ShapeFunctionPool.cpp

@@ -0,0 +1,172 @@
+/*! \file ShapeFunctionPool.cpp
+    \brief Compute shape functions and their gradients with respect to
+     the local coodinates, and store the results
+
+     \author Wenqing Wang
+     \date Feb. 2015
+
+     \copyright
+      Copyright (c) 2015, OpenGeoSys Community (http://www.opengeosys.org)
+             Distributed under a Modified BSD License.
+             See accompanying file LICENSE.txt or
+             http://www.opengeosys.org/project/license
+*/
+
+#include "ShapeFunctionPool.h"
+
+#include <cassert> /* assert */
+#include "fem_ele.h"
+
+namespace FiniteElement
+{
+ShapeFunctionPool::ShapeFunctionPool(const std::vector<MshElemType::type>& elem_types, CElement& quadrature,
+                                     const int num_sample_gs_pnts)
+{
+	int num_elem_nodes[2][MshElemType::NUM_ELEM_TYPES];
+	int dim_elem[MshElemType::NUM_ELEM_TYPES];
+
+	int id = static_cast<int>(MshElemType::LINE) - 1;
+	num_elem_nodes[0][id] = 2;
+	num_elem_nodes[1][id] = 3;
+	dim_elem[id] = 1;
+
+	id = static_cast<int>(MshElemType::QUAD) - 1;
+	num_elem_nodes[0][id] = 4;
+	num_elem_nodes[1][id] = 9;
+	dim_elem[id] = 2;
+
+	id = static_cast<int>(MshElemType::QUAD8) - 1;
+	num_elem_nodes[0][id] = 4;
+	num_elem_nodes[1][id] = 8;
+	dim_elem[id] = 2;
+
+	id = static_cast<int>(MshElemType::TRIANGLE) - 1;
+	num_elem_nodes[0][id] = 3;
+	num_elem_nodes[1][id] = 6;
+	dim_elem[id] = 2;
+
+	id = static_cast<int>(MshElemType::HEXAHEDRON) - 1;
+	num_elem_nodes[0][id] = 8;
+	num_elem_nodes[1][id] = 20;
+	dim_elem[id] = 3;
+
+	id = static_cast<int>(MshElemType::TETRAHEDRON) - 1;
+	num_elem_nodes[0][id] = 4;
+	num_elem_nodes[1][id] = 10;
+	dim_elem[id] = 3;
+
+	id = static_cast<int>(MshElemType::PRISM) - 1;
+	num_elem_nodes[0][id] = 6;
+	num_elem_nodes[1][id] = 15;
+	dim_elem[id] = 3;
+
+	id = static_cast<int>(MshElemType::PYRAMID) - 1;
+	num_elem_nodes[0][id] = 5;
+	num_elem_nodes[1][id] = 13;
+	dim_elem[id] = 3;
+
+	// std::vector<int> elem_type_ids
+	const std::size_t n_ele_types = elem_types.size();
+	_shape_function.resize(n_ele_types);
+	_shape_function_center.resize(n_ele_types);
+	_grad_shape_function.resize(n_ele_types);
+	_grad_shape_function_center.resize(n_ele_types);
+	for (std::size_t i = 0; i < elem_types.size(); i++)
+	{
+		const MshElemType::type e_type = elem_types[i];
+		if (e_type == MshElemType::INVALID)
+			continue;
+		// Set number of integration points.
+		quadrature.SetGaussPointNumber(num_sample_gs_pnts);
+		quadrature.SetIntegrationPointNumber(e_type);
+
+		const int type_id = static_cast<int>(e_type) - 1;
+		int num_int_pnts = quadrature.GetNumGaussPoints();
+		const int num_nodes = num_elem_nodes[quadrature.getOrder() - 1][type_id];
+
+		std::vector<double> elem_shape_function_center(num_nodes);
+		_shape_function_center[type_id] = elem_shape_function_center;
+
+		const int size_shape_fct = num_nodes * num_int_pnts;
+		std::vector<double> elem_shape_function(size_shape_fct);
+		_shape_function[type_id] = elem_shape_function;
+
+		std::vector<double> elem_grad_shape_function(dim_elem[type_id] * size_shape_fct);
+		_grad_shape_function[type_id] = elem_grad_shape_function;
+
+		std::vector<double> elem_grad_shape_function_center(dim_elem[type_id] * num_nodes);
+		_grad_shape_function_center[type_id] = elem_grad_shape_function_center;
+	}
+
+	computeQuadratures(elem_types, num_elem_nodes, dim_elem, quadrature, num_sample_gs_pnts);
+}
+
+void ShapeFunctionPool::computeQuadratures(const std::vector<MshElemType::type>& elem_types,
+                                           const int num_elem_nodes[2][MshElemType::NUM_ELEM_TYPES],
+                                           const int dim_elem[],
+                                           CElement& quadrature, const int num_sample_gs_pnts)
+{
+	const int order = quadrature.getOrder();
+	for (std::size_t i = 0; i < elem_types.size(); i++)
+	{
+		const MshElemType::type e_type = elem_types[i];
+		if (e_type == MshElemType::INVALID)
+			continue;
+
+		const int type_id = static_cast<int>(e_type) - 1;
+		quadrature.ConfigShapefunction(e_type);
+
+		const int nnodes = num_elem_nodes[order - 1][type_id];
+		const int elem_dim = dim_elem[type_id];
+
+		double* shape_function_center_values = _shape_function_center[type_id].data();
+		quadrature.SetCenterGP(e_type);
+		quadrature.ComputeShapefct(order, shape_function_center_values);
+		double* grad_shape_function_center_values = _grad_shape_function_center[type_id].data();
+		quadrature.computeGradShapefctLocal(order, grad_shape_function_center_values);
+
+		double* shape_function_values = _shape_function[type_id].data();
+		double* dshape_function_values = _grad_shape_function[type_id].data();
+		// Set number of integration points.
+		quadrature.SetGaussPointNumber(num_sample_gs_pnts);
+		quadrature.SetIntegrationPointNumber(e_type);
+
+		for (int gp = 0; gp < quadrature.GetNumGaussPoints(); gp++)
+		{
+			int gp_r, gp_s, gp_t;
+			quadrature.SetGaussPoint(e_type, gp, gp_r, gp_s, gp_t);
+			double* shape_function_values_gs = &shape_function_values[gp * nnodes];
+			quadrature.ComputeShapefct(order, shape_function_values_gs);
+
+			double* dshape_function_values_gs = &dshape_function_values[gp * nnodes * elem_dim];
+			quadrature.computeGradShapefctLocal(order, dshape_function_values_gs);
+		}
+	}
+}
+
+const double* ShapeFunctionPool::getShapeFunctionValues(const MshElemType::type elem_type) const
+{
+	return _shape_function[static_cast<int>(elem_type) - 1].data();
+}
+
+unsigned ShapeFunctionPool::getShapeFunctionArraySize(const MshElemType::type elem_type) const
+{
+	return _shape_function[static_cast<int>(elem_type) - 1].size();
+}
+
+const double* ShapeFunctionPool::getShapeFunctionCenterValues(const MshElemType::type elem_type) const
+{
+	return _shape_function_center[static_cast<int>(elem_type) - 1].data();
+}
+
+const double* ShapeFunctionPool::getGradShapeFunctionValues(const MshElemType::type elem_type) const
+{
+	return _grad_shape_function[static_cast<int>(elem_type) - 1].data();
+}
+
+const double* ShapeFunctionPool::getGradShapeFunctionCenterValues(const MshElemType::type elem_type) const
+{
+	return _grad_shape_function_center[static_cast<int>(elem_type) - 1].data();
+}
+
+} // end namespace

FEM/ShapeFunctionPool.h

@@ -0,0 +1,73 @@
+/*! \file ShapeFunctionPool.h
+    \brief Compute shape functions and their gradients with respect to
+     the local coodinates, and store the results
+
+     \author Wenqing Wang
+     \date Feb. 2015
+
+     \copyright
+      Copyright (c) 2015, OpenGeoSys Community (http://www.opengeosys.org)
+             Distributed under a Modified BSD License.
+             See accompanying file LICENSE.txt or
+             http://www.opengeosys.org/project/license
+*/
+#ifndef OGS_SHAPEFUNCTIONPOOL_H
+#define OGS_SHAPEFUNCTIONPOOL_H
+
+#include <vector>
+
+#include "MSHEnums.h"
+
+namespace FiniteElement
+{
+class CElement;
+
+class ShapeFunctionPool
+{
+public:
+	/*!
+	    \param elem_types          All involved element types.
+	    \param quadrature          Numerical integration object.
+	    \param num_sample_gs_pnts  Number of sample Gauss points.
+	*/
+	ShapeFunctionPool(const std::vector<MshElemType::type>& elem_types, CElement& quadrature,
+	                  const int num_sample_gs_pnts);
+	~ShapeFunctionPool() {}
+
+	/// Get shape function values of an element type
+	const double* getShapeFunctionValues(const MshElemType::type elem_type) const;
+	/// Get the size of shape function array of an element type.
+	unsigned getShapeFunctionArraySize(const MshElemType::type elem_type) const;
+
+	/// Get shape function values at the element centroid of an element type
+	const double* getShapeFunctionCenterValues(const MshElemType::type elem_type) const;
+
+	/// Get the values of the gradient of shape function of an element type
+	const double* getGradShapeFunctionValues(const MshElemType::type elem_type) const;
+
+	/// Get the gradient of shape function values at the element centroid of an element type
+	const double* getGradShapeFunctionCenterValues(const MshElemType::type elem_type) const;
+
+private:
+	/// Results of shape functions of all integration points.
+	std::vector< std::vector<double> > _shape_function;
+
+	/// Results of shape functions of all integration points at element centroid.
+	std::vector< std::vector<double> > _shape_function_center;
+
+	/// Results of the gradient of shape functions with respect to
+	/// local coordinates of all integration points.
+	std::vector< std::vector<double> > _grad_shape_function;
+
+	/// Results of the gradient of shape functions of all integration points at element centroid.
+	std::vector< std::vector<double> > _grad_shape_function_center;
+
+	void computeQuadratures(const std::vector<MshElemType::type>& elem_types,
+	                        const int num_elem_nodes[2][MshElemType::NUM_ELEM_TYPES],
+	                        const int dim_elem[],
+	                        CElement& quadrature,
+	                        const int num_sample_gs_pnts);
+};
+} // end namespace
+
+#endif