Skip to content

Commit

Permalink
Merge pull request #38 from wenqing/shapefct_cache
Browse files Browse the repository at this point in the history
Improvement in shape function computation
  • Loading branch information
norihiro-w authored Jun 24, 2016
2 parents 407a0f5 + c1feb7d commit 730a841
Show file tree
Hide file tree
Showing 37 changed files with 1,987 additions and 1,226 deletions.
2 changes: 2 additions & 0 deletions FEM/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ set( HEADERS
tools.h
vtk.h
OutputTools.h
ShapeFunctionPool.h
)

set( SOURCES
Expand Down Expand Up @@ -102,6 +103,7 @@ set( SOURCES
Stiff_Bulirsch-Stoer.cpp
tools.cpp
vtk.cpp
ShapeFunctionPool.cpp
)

if(NOT OGS_LSOLVER STREQUAL PETSC)
Expand Down
2 changes: 1 addition & 1 deletion FEM/DUMUX.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 = "";
Expand Down
23 changes: 12 additions & 11 deletions FEM/Output.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include <string>

#include <cfloat> // DBL_EPSILON
#include <algorithm> // remove-if

#include "BuildInfo.h"
#include "FEMIO/GeoIO.h"
Expand Down Expand Up @@ -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);
}

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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++)
{
Expand All @@ -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]);
Expand All @@ -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;
}

Expand Down
6 changes: 2 additions & 4 deletions FEM/OutputTools.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
172 changes: 172 additions & 0 deletions FEM/ShapeFunctionPool.cpp
Original file line number Diff line number Diff line change
@@ -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
73 changes: 73 additions & 0 deletions FEM/ShapeFunctionPool.h
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 730a841

Please sign in to comment.