v0.14.0
Loading...
Searching...
No Matches
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.
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
struct OpVolume : public VolumeElementForcesAndSourcesCore::UserDataOperator {
Vec vOl;
OpVolume(const std::string &field_name, Vec vol)
: VolumeElementForcesAndSourcesCore::UserDataOperator(field_name, OPROW),
vOl(vol) {}
if (type != MBVERTEX)
const int nb_int_pts = getGaussPts().size2();
// cerr << nb_int_pts << endl;
auto t_w = getFTensor0IntegrationWeight();
double v = getMeasure();
double vol = 0;
for (int gg = 0; gg != nb_int_pts; gg++) {
vol += t_w * v;
++t_w;
}
CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
// PetscPrintf(PETSC_COMM_WORLD,"domain: calculate matrix\n");
}
};
struct OpFace : public FaceElementForcesAndSourcesCore::UserDataOperator {
Vec vOl;
OpFace(const std::string &field_name, Vec vol)
vOl(vol) {}
if (type != MBVERTEX)
const int nb_int_pts = getGaussPts().size2();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_coords = getFTensor1CoordsAtGaussPts();
FTensor::Index<'i', 3> i;
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;
CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
}
MoFEMErrorCode doWork(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;
CHKERR m_field.getInterface(simple_interface);
{
// get options from command line
CHKERR simple_interface->getOptions();
// load mesh file
CHKERR simple_interface->loadFile();
// add fields
CHKERR simple_interface->addDomainField("MESH_NODE_POSITIONS", H1,
CHKERR simple_interface->addBoundaryField("MESH_NODE_POSITIONS", H1,
CHKERR simple_interface->addSkeletonField("MESH_NODE_POSITIONS", H1,
// set fields order
CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS", 1);
// setup problem
CHKERR simple_interface->setUp();
// Project mesh coordinate on mesh
Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
// get dm
auto dm = simple_interface->getDM();
// create elements
auto domain_fe =
boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
auto boundary_fe =
boost::make_shared<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};
auto vol = createGhostVector(
PETSC_COMM_WORLD, m_field.get_comm_rank() == 0 ? 1 : 0, 1,
m_field.get_comm_rank() == 0 ? 0 : 1, ghosts);
auto surf_vol = vectorDuplicate(vol);
// set operator to the volume element
auto material_grad_mat = boost::make_shared<MatrixDouble>();
auto material_det_vec = boost::make_shared<VectorDouble>();
domain_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
material_grad_mat));
domain_fe->getOpPtrVector().push_back(new OpInvertMatrix<3>(
material_grad_mat, material_det_vec, nullptr));
domain_fe->getOpPtrVector().push_back(
new OpSetHOWeights(material_det_vec));
domain_fe->getOpPtrVector().push_back(
new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
domain_fe->getOpPtrVector().push_back(
new OpVolume("MESH_NODE_POSITIONS", vol));
// set operator to the face element
boundary_fe->getOpPtrVector().push_back(
new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
boundary_fe->getOpPtrVector().push_back(
new OpFace("MESH_NODE_POSITIONS", surf_vol));
// make integration in volume (here real calculations starts)
domain_fe);
// make integration on boundary
boundary_fe);
auto skeleton_fe = boost::make_shared<FEMethod>();
auto A = createDMMatrix(dm);
auto B = createDMMatrix(dm);
auto f = createDMVector(dm);
auto x = createDMVector(dm);
auto x_t = createDMVector(dm);
auto x_tt = createDMVector(dm);
skeleton_fe->f = f;
skeleton_fe->A = A;
skeleton_fe->B = B;
skeleton_fe->x = x;
skeleton_fe->x_t = x_t;
skeleton_fe->x_tt = x_tt;
skeleton_fe->preProcessHook = [&]() {
if (f != skeleton_fe->ts_F)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
if (A != skeleton_fe->ts_A)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
if (B != skeleton_fe->ts_B)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
if (x != skeleton_fe->ts_u)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
if (x_t != skeleton_fe->ts_u_t)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
if (x_tt != skeleton_fe->ts_u_tt)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
};
skeleton_fe->postProcessHook = []() { return 0; };
skeleton_fe->operatorHook = []() { return 0; };
skeleton_fe);
// assemble volumes from processors and accumulate on processor of rank 0
CHKERR VecAssemblyBegin(vol);
CHKERR VecAssemblyEnd(vol);
CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(surf_vol);
CHKERR VecAssemblyEnd(surf_vol);
CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
if (m_field.get_comm_rank() == 0) {
double *a_vol;
CHKERR VecGetArray(vol, &a_vol);
double *a_surf_vol;
CHKERR VecGetArray(surf_vol, &a_surf_vol);
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");
}
CHKERR VecRestoreArray(vol, &a_vol);
CHKERR VecRestoreArray(vol, &a_surf_vol);
}
}
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static char help[]
int main()
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:988
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
constexpr AssemblyType A
constexpr auto field_name
Set integration rule to boundary elements.
int operator()(int, int, int) const
virtual int get_comm_rank() const =0
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Calculate HO coordinates at gauss points.
Set inverse jacobian to base functions.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition Simple.hpp:348
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition Simple.hpp:355
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:264
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:473
MoFEMErrorCode addBoundaryField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition Simple.cpp:282
MoFEMErrorCode addSkeletonField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on skeleton.
Definition Simple.cpp:300
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:611
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:341
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
OpFace(const Range &skin_ents, const std::vector< unsigned char > &marker)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpVolume(const std::string &field_name, Vec vol)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Set integration rule.
int operator()(int, int, int) const