v0.6.10
simple_interface.cpp

Calculate volume by integrating volume elements and using divergence theorem by integrating surface elements.

/**
* \file simple_interface.cpp
* \ingroup mofem_simple_interface
* \example simple_interface.cpp
*
* Calculate volume by integrating volume elements and using divergence theorem
* by integrating surface elements.
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
Vec vOl;
OpVolume(const std::string &field_name,Vec vol):
vOl(vol) {
}
MoFEMErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
if(type!=MBVERTEX) MoFEMFunctionReturnHot(0);
const int nb_int_pts = getGaussPts().size2();
// cerr << nb_int_pts << endl;
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
FTensor::Tensor0<double*> t_ho_det = getFTenosr0HoMeasure();
double v = getMeasure();
double vol = 0;
for(int gg = 0;gg!=nb_int_pts;gg++) {
vol += t_w*t_ho_det*v;
// cerr << t_ho_det << endl;
++t_w;
++t_ho_det;
}
ierr = VecSetValue(vOl,0,vol,ADD_VALUES); CHKERRG(ierr);
}
int row_side,int col_side,
EntityType row_type, EntityType col_type,
) {
// PetscPrintf(PETSC_COMM_WORLD,"domain: calculate matrix\n");
}
};
Vec vOl;
OpFace(const std::string &field_name,Vec vol):
vOl(vol) {
}
MoFEMErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
if(type!=MBVERTEX) MoFEMFunctionReturnHot(0);
const int nb_int_pts = getGaussPts().size2();
FTensor::Tensor1<double*,3> t_normal = getTensor1NormalsAtGaussPt();
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
FTensor::Tensor1<double*,3> t_coords = getTensor1HoCoordsAtGaussPts();
double vol = 0;
for(int gg = 0;gg!=nb_int_pts;gg++) {
vol += (t_coords(i)*t_normal(i))*t_w;
++t_normal;
++t_w;
++t_coords;
}
vol /= 6;
ierr = VecSetValue(vOl,0,vol,ADD_VALUES); CHKERRG(ierr);
}
int row_side,int col_side,
EntityType row_type, EntityType col_type,
) {
// PetscPrintf(PETSC_COMM_WORLD,"boundary: calculate matrix\n");
}
};
struct VolRule { int operator()(int,int,int) const { return 2; } };
struct FaceRule { int operator()(int,int,int) const { return 4; } };
int main(int argc, char *argv[]) {
//
// initialize petsc
MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
try {
// Create MoAB database
moab::Core moab_core;
moab::Interface& moab = moab_core;
// Create MoFEM database and link it to MoAB
MoFEM::Core mofem_core(moab);
MoFEM::Interface& m_field = mofem_core;
// Register DM Manager
DMType dm_name = "DMMOFEM";
// Simple interface
Simple *simple_interface;
ierr = m_field.getInterface(simple_interface); CHKERRG(ierr);
{
// get options from command line
ierr = simple_interface->getOptions(); CHKERRG(ierr);
// load mesh file
ierr = simple_interface->loadFile(); CHKERRG(ierr);
// add fields
ierr = simple_interface->addDomainField("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRG(ierr);
ierr = simple_interface->addBoundaryField("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRG(ierr);
// set fields order
ierr = simple_interface->setFieldOrder("MESH_NODE_POSITIONS",1); CHKERRG(ierr);
// setup problem
ierr = simple_interface->setUp(); CHKERRG(ierr);
// Project mesh coordinate on mesh
Projection10NodeCoordsOnField ent_method(m_field,"MESH_NODE_POSITIONS");
ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method); CHKERRG(ierr);
DM dm;
// get dm
ierr = simple_interface->getDM(&dm); CHKERRG(ierr);
// create elements
boost::shared_ptr<ForcesAndSourcesCore> domain_fe =
boost::shared_ptr<ForcesAndSourcesCore>(new VolumeElementForcesAndSourcesCore(m_field));
boost::shared_ptr<ForcesAndSourcesCore> boundary_fe =
boost::shared_ptr<ForcesAndSourcesCore>(new FaceElementForcesAndSourcesCore(m_field));
// set integration rule
domain_fe->getRuleHook = VolRule();
boundary_fe->getRuleHook = FaceRule();
// create distributed vector to accumulate values from processors.
int ghosts[] = { 0 };
Vec vol,surf_vol;
ierr = VecCreateGhost(
PETSC_COMM_WORLD,m_field.get_comm_rank()==0?1:0,1,m_field.get_comm_rank()==0?0:1,ghosts,&vol
ierr = VecDuplicate(vol,&surf_vol); CHKERRG(ierr);
// set operator to the volume element
domain_fe->getOpPtrVector().push_back(new OpVolume("MESH_NODE_POSITIONS",vol)
);
// set operator to the face element
boundary_fe->getOpPtrVector().push_back(new OpFace("MESH_NODE_POSITIONS",surf_vol)
);
// make integration in volume (here real calculations starts)
ierr = DMoFEMLoopFiniteElements(dm,simple_interface->getDomainFEName(),domain_fe); CHKERRG(ierr);
// make integration on boundary
ierr = DMoFEMLoopFiniteElements(dm,simple_interface->getBoundaryFEName(),boundary_fe); CHKERRG(ierr);
// assemble volumes from processors and accumulate on processor of rank 0
ierr = VecAssemblyBegin(vol); CHKERRG(ierr);
ierr = VecAssemblyEnd(vol); CHKERRG(ierr);
ierr = VecGhostUpdateBegin(vol,ADD_VALUES,SCATTER_REVERSE); CHKERRG(ierr);
ierr = VecGhostUpdateEnd(vol,ADD_VALUES,SCATTER_REVERSE); CHKERRG(ierr);
ierr = VecAssemblyBegin(surf_vol); CHKERRG(ierr);
ierr = VecAssemblyEnd(surf_vol); CHKERRG(ierr);
ierr = VecGhostUpdateBegin(surf_vol,ADD_VALUES,SCATTER_REVERSE); CHKERRG(ierr);
ierr = VecGhostUpdateEnd(surf_vol,ADD_VALUES,SCATTER_REVERSE); CHKERRG(ierr);
if(m_field.get_comm_rank()==0) {
double *a_vol;
ierr = VecGetArray(vol,&a_vol); CHKERRG(ierr);
double *a_surf_vol;
ierr = VecGetArray(surf_vol,&a_surf_vol); CHKERRG(ierr);
cout << "Volume = " << a_vol[0] << endl;
cout << "Surf Volume = " << a_surf_vol[0] << endl;
if(fabs(a_vol[0]-a_surf_vol[0])>1e-12) {
SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"Should be zero");
}
ierr = VecRestoreArray(vol,&a_vol); CHKERRG(ierr);
ierr = VecRestoreArray(vol,&a_surf_vol); CHKERRG(ierr);
}
// destroy vector
ierr = VecDestroy(&vol); CHKERRG(ierr);
ierr = VecDestroy(&surf_vol); CHKERRG(ierr);
// destroy dm
ierr = DMDestroy(&dm); CHKERRG(ierr);
}
} catch (MoFEMException const &e) {
SETERRQ(PETSC_COMM_SELF,e.errorCode,e.errorMessage);
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}