v0.15.0
Loading...
Searching...
No Matches
mofem/tutorials/scl-12/electrostatics.cpp
/**
* @file Electrostatics.cpp
* \example mofem/tutorials/scl-12/electrostatics.cpp
* */
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
static char help[] = "...\n\n";
public:
// Declaration of the main function to run analysis
private:
// Declaration of other main functions called in runProgram()
int oRder = 2; // default order
int geomOrder = 1; // default gemoetric order
Simple *simpleInterface;
std::string domainField;
boost::shared_ptr<std::map<int, BlockData>> permBlockSetsPtr;
boost::shared_ptr<std::map<int, BlockData>> intBlockSetsPtr;
boost::shared_ptr<std::map<int, BlockData>> electrodeBlockSetsPtr;
boost::shared_ptr<DataAtIntegrationPts> commonDataPtr;
boost::shared_ptr<ForcesAndSourcesCore> interFaceRhsFe;
boost::shared_ptr<ForcesAndSourcesCore> electrodeRhsFe;
double aLpha = 0.0; // declaration for total charge on first electrode
double bEta = 0.0; // declaration for total charge on second electrode
SmartPetscObj<Vec> petscVec; // petsc vector for the charge solution
SmartPetscObj<Vec> petscVecEnergy; // petsc vector for the energy solution
PetscBool out_skin = PETSC_FALSE; //
PetscBool is_partitioned = PETSC_FALSE;
enum VecElements { ZERO = 0, ONE = 1, LAST_ELEMENT };
int atomTest = 0;
};
: domainField("POTENTIAL"), mField(m_field) {}
//! [Read mesh]
true; // create lower dimensional element to ensure sharing of partitioed
// entities at the boundaries.
}
//! [Read mesh]
//! [Setup problem]
Range domain_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
true);
auto get_ents_by_dim = [&](const auto dim) {
if (dim == SPACE_DIM) {
return domain_ents;
} else {
Range ents;
if (dim == 0)
CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
else
CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
return ents;
}
};
// Select base for the field based on the element type
auto get_base = [&]() {
auto domain_ents = get_ents_by_dim(SPACE_DIM);
if (domain_ents.empty())
const auto type = type_from_handle(domain_ents[0]);
switch (type) {
case MBQUAD:
case MBHEX:
case MBTRI:
case MBTET:
default:
CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type is not handled");
}
return NOBASE;
};
auto base = get_base();
CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &oRder, PETSC_NULLPTR);
CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geomOrder,
PETSC_NULLPTR);
CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atomTest,
PETSC_NULLPTR);
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
return mField.loop_dofs("GEOMETRY", ent_method);
};
CHKERR project_ho_geometry();
commonDataPtr = boost::make_shared<DataAtIntegrationPts>(mField);
// gets the map of the permittivity attributes and the block sets
permBlockSetsPtr = boost::make_shared<std::map<int, BlockData>>();
Range mat_electr_ents; // range of entities with the permittivity
if (bit->getName().compare(0, 12, "MAT_ELECTRIC") == 0) {
const int id = bit->getMeshsetId();
auto &block_data = (*permBlockSetsPtr)[id];
CHKERR mField.get_moab().get_entities_by_dimension(
bit->getMeshset(), SPACE_DIM, block_data.domainEnts, true);
mat_electr_ents.merge(block_data.domainEnts);
std::vector<double> attributes;
bit->getAttributes(attributes);
if (attributes.size() < 1) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
" At least one permittivity attributes should be given but "
"found %d",
attributes.size());
}
block_data.iD = id; // id of the block
block_data.epsPermit = attributes[0]; // permittivity value of the block
}
}
// gets the map of the charge attributes and the block sets
intBlockSetsPtr = boost::make_shared<std::map<int, BlockData>>();
Range int_electr_ents; // range of entities with the charge
if (bit->getName().compare(0, 12, "INT_ELECTRIC") == 0) {
const int id = bit->getMeshsetId();
auto &block_data = (*intBlockSetsPtr)[id];
CHKERR mField.get_moab().get_entities_by_dimension(
bit->getMeshset(), SPACE_DIM - 1, block_data.interfaceEnts, true);
int_electr_ents.merge(block_data.interfaceEnts);
std::vector<double> attributes;
bit->getAttributes(attributes);
if (attributes.size() < 1) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"At least one charge attributes should be given but found %d",
attributes.size());
}
block_data.iD = id; // id-> block ID
block_data.chargeDensity = attributes[0]; // block charge attribute
}
}
// gets the map of the electrode entity range in the block sets
electrodeBlockSetsPtr = boost::make_shared<std::map<int, BlockData>>();
Range electrode_ents; // range of entities with the electrode
int electrodeCount = 0;
if (bit->getName().compare(0, 9, "ELECTRODE") == 0) {
const int id = bit->getMeshsetId();
auto &block_data = (*electrodeBlockSetsPtr)[id];
++electrodeCount;
CHKERR mField.get_moab().get_entities_by_dimension(
bit->getMeshset(), SPACE_DIM - 1, block_data.electrodeEnts, true);
electrode_ents.merge(block_data.electrodeEnts);
if (electrodeCount > 2) {
SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"Three or more electrode blocksets found");
;
}
}
}
// sync entities
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
mat_electr_ents);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
int_electr_ents);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
electrode_ents);
CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-out_skin", &out_skin,
PETSC_NULLPTR);
CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_partitioned", &is_partitioned,
PETSC_NULLPTR);
// get the skin entities
Skinner skinner(&mField.get_moab());
Range skin_tris;
CHKERR skinner.find_skin(0, mat_electr_ents, false, skin_tris);
Range proc_skin;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
CHKERR pcomm->filter_pstatus(skin_tris,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &proc_skin);
} else {
proc_skin = skin_tris;
}
// add the skin entities to the field
"SKIN");
// add the interface entities to the field
SPACE_DIM - 1, "INTERFACE");
// add the electrode entities to the field
"ELECTRODE");
// sync field entities
mField.getInterface<CommInterface>()->synchroniseFieldEntities(domainField);
mField.getInterface<CommInterface>()->synchroniseFieldEntities("GEOMETRY");
DMType dm_name = "DMMOFEM";
SmartPetscObj<DM> dm;
dm = createDM(mField.get_comm(), dm_name);
// create dm instance
CHKERR DMSetType(dm, dm_name);
// initialise petsc vector for required processor
int local_size;
if (mField.get_comm_rank() == 0) // get_comm_rank() gets processor number
local_size = LAST_ELEMENT; // last element gives size of vector
else
// other processors (e.g. 1, 2, 3, etc.)
local_size = 0; // local size of vector is zero on other processors
}
//! [Setup problem]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
// Remove_BCs_from_blockset name "BOUNDARY_CONDITION";
CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
simpleInterface->getProblemName(), "BOUNDARY_CONDITION",
std::string(domainField), true);
}
//! [Boundary condition]
//! [Set integration rules]
auto rule_lhs = [this](int, int, int p) -> int { return 2 * p + geomOrder -1; };
auto rule_rhs = [this](int, int, int p) -> int { return 2 * p + geomOrder -1; };
auto pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
}
//! [Set integration rules]
//! [Assemble system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
commonDataPtr = boost::make_shared<DataAtIntegrationPts>(mField);
auto add_domain_base_ops = [&](auto &pipeline) {
"GEOMETRY");
pipeline.push_back(
};
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
auto epsilon = [&](const double, const double, const double) {
return commonDataPtr->blockPermittivity;
};
{ // Push operators to the Pipeline that is responsible for calculating LHS
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating RHS
auto set_values_to_bc_dofs = [&](auto &fe) {
auto get_bc_hook = [&]() {
EssentialPreProc<TemperatureCubitBcData> hook(mField, fe, {});
return hook;
};
fe->preProcessHook = get_bc_hook();
};
// Set essential BC
auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
using OpInternal =
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>(domainField,
grad_u_ptr));
auto minus_epsilon = [&](double, double, double) constexpr {
return -commonDataPtr->blockPermittivity;
};
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInternal(domainField, grad_u_ptr, minus_epsilon));
};
set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
calculate_residual_from_set_values_on_bc(
pipeline_mng->getOpDomainRhsPipeline());
auto bodySourceTerm = [&](const double, const double, const double) {
return bodySource;
};
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBodySourceVectorb(domainField, bodySourceTerm));
}
interFaceRhsFe = boost::shared_ptr<ForcesAndSourcesCore>(
interFaceRhsFe->getRuleHook = [this](int, int, int p) {
return 2 * p + geomOrder -1;
};
{
interFaceRhsFe->getOpPtrVector(), {NOSPACE}, "GEOMETRY");
interFaceRhsFe->getOpPtrVector().push_back(
auto sIgma = [&](const double, const double, const double) {
return commonDataPtr->blockChrgDens;
};
interFaceRhsFe->getOpPtrVector().push_back(
}
}
//! [Assemble system]
//! [Solve system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto ksp_solver = pipeline_mng->createKSP();
boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element does
DM dm;
CHKERR DMMoFEMKSPSetComputeRHS(dm, "INTERFACE", interFaceRhsFe, null, null);
CHKERR KSPSetFromOptions(ksp_solver);
// Create RHS and solution vectors
auto F = createDMVector(dm);
auto D = vectorDuplicate(F);
// Solve the system
CHKERR KSPSetUp(ksp_solver);
CHKERR KSPSolve(ksp_solver, F, D);
CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
double fnorm;
CHKERR VecNorm(F, NORM_2, &fnorm);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "F norm = %9.8e\n", fnorm);
// Scatter result data on the mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
double dnorm;
CHKERR VecNorm(D, NORM_2, &dnorm);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "D norm = %9.8e\n", dnorm);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve system]
//! [Output results]
auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset(); //
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
// lamda function to calculate electric field
auto calculate_e_field = [&](auto &pipeline) {
auto u_ptr = boost::make_shared<VectorDouble>();
auto x_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto e_field_ptr = boost::make_shared<MatrixDouble>();
// add higher order operator
"GEOMETRY");
// calculate field values
pipeline.push_back(new OpCalculateScalarFieldValues(domainField, u_ptr));
pipeline.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
// calculate gradient
pipeline.push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>(domainField, grad_u_ptr));
// calculate electric field
pipeline.push_back(new OpElectricField(e_field_ptr, grad_u_ptr));
return boost::make_tuple(u_ptr, e_field_ptr, x_ptr);
};
auto [u_ptr, e_field_ptr, x_ptr] =
calculate_e_field(post_proc_fe->getOpPtrVector());
auto e_field_times_perm_ptr = boost::make_shared<MatrixDouble>();
auto energy_density_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpGradTimesPerm(domainField, e_field_ptr, e_field_times_perm_ptr,
post_proc_fe->getOpPtrVector().push_back(
new OpEnergyDensity(domainField, e_field_ptr, energy_density_ptr,
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"POTENTIAL", u_ptr},
{"ENERGY_DENSITY", energy_density_ptr}},
OpPPMap::DataMapMat{
{"GEOMETRY", x_ptr},
{"ELECTRIC_FIELD", e_field_ptr},
{"ELECTRIC_DISPLACEMENT", e_field_times_perm_ptr},
},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{})
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out.h5m");
if (out_skin && SPACE_DIM == 3) {
auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
auto op_loop_skin = new OpLoopSide<SideEle>(
mField, simpleInterface->getDomainFEName(), SPACE_DIM);
auto [u_ptr, e_field_ptr, x_ptr] =
calculate_e_field(op_loop_skin->getOpPtrVector());
op_loop_skin->getOpPtrVector().push_back(
new OpGradTimesPerm(domainField, e_field_ptr, e_field_times_perm_ptr,
permBlockSetsPtr, commonDataPtr));
op_loop_skin->getOpPtrVector().push_back(
new OpEnergyDensity(domainField, e_field_ptr, energy_density_ptr,
permBlockSetsPtr, commonDataPtr));
// push op to boundary element
post_proc_skin->getOpPtrVector().push_back(op_loop_skin);
post_proc_skin->getOpPtrVector().push_back(new OpPPMap(
post_proc_skin->getPostProcMesh(), post_proc_skin->getMapGaussPts(),
OpPPMap::DataMapVec{{"POTENTIAL", u_ptr},
{"ENERGY_DENSITY", energy_density_ptr}},
OpPPMap::DataMapMat{{"ELECTRIC_FIELD", e_field_ptr},
{"GEOMETRY", x_ptr},
{"ELECTRIC_DISPLACEMENT", e_field_times_perm_ptr}},
OpPPMap::DataMapMat{}, OpPPMap::DataMapMat{}));
CHKERR DMoFEMLoopFiniteElements(simpleInterface->getDM(), "SKIN",
post_proc_skin);
CHKERR post_proc_skin->writeFile("out_skin.h5m");
}
}
//! [Output results]
//! [Get Total Energy]
auto pip_energy = mField.getInterface<PipelineManager>();
pip_energy->getDomainRhsFE().reset();
pip_energy->getDomainLhsFE().reset();
pip_energy->getBoundaryLhsFE().reset();
pip_energy->getBoundaryRhsFE().reset();
pip_energy->getOpDomainRhsPipeline().clear();
pip_energy->getOpDomainLhsPipeline().clear();
// gets the map of the internal domain entity range to get the total energy
boost::shared_ptr<std::map<int, BlockData>> intrnlDomnBlckSetPtr =
boost::make_shared<std::map<int, BlockData>>();
Range internal_domain; // range of entities marked the internal domain
if (bit->getName().compare(0, 10, "DOMAIN_INT") == 0) {
const int id = bit->getMeshsetId();
auto &block_data = (*intrnlDomnBlckSetPtr)[id];
CHKERR mField.get_moab().get_entities_by_dimension(
bit->getMeshset(), SPACE_DIM, block_data.internalDomainEnts, true);
internal_domain.merge(block_data.internalDomainEnts);
block_data.iD = id;
}
}
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
internal_domain);
pip_energy->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto e_field_ptr = boost::make_shared<MatrixDouble>();
pip_energy->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>(domainField, grad_u_ptr));
pip_energy->getOpDomainRhsPipeline().push_back(
new OpElectricField(e_field_ptr, grad_u_ptr));
commonDataPtr = boost::make_shared<DataAtIntegrationPts>(mField);
pip_energy->getOpDomainRhsPipeline().push_back(
intrnlDomnBlckSetPtr, commonDataPtr, petscVecEnergy));
pip_energy->loopFiniteElements();
CHKERR VecAssemblyBegin(petscVecEnergy);
CHKERR VecAssemblyEnd(petscVecEnergy);
double total_energy = 0.0; // declaration for total energy
const double *array;
CHKERR VecGetArrayRead(petscVecEnergy, &array);
total_energy = array[ZERO];
MOFEM_LOG_C("SELF", Sev::inform, "Total Energy: %6.15f", total_energy);
CHKERR VecRestoreArrayRead(petscVecEnergy, &array);
}
}
//! [Get Total Energy]
//! [Get Charges]
auto op_loop_side = new OpLoopSide<SideEle>(
op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
auto grad_u_ptr_charge = boost::make_shared<MatrixDouble>();
auto e_ptr_charge = boost::make_shared<MatrixDouble>();
op_loop_side->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>(domainField,
grad_u_ptr_charge));
op_loop_side->getOpPtrVector().push_back(
new OpElectricField(e_ptr_charge, grad_u_ptr_charge));
auto d_jump = boost::make_shared<MatrixDouble>();
op_loop_side->getOpPtrVector().push_back(new OpElectricDispJump<SPACE_DIM>(
domainField, e_ptr_charge, d_jump, commonDataPtr, permBlockSetsPtr));
electrodeRhsFe = boost::shared_ptr<ForcesAndSourcesCore>(
electrodeRhsFe->getRuleHook = [this](int, int, int p) {
return 2 * p + geomOrder -1;
};
// push all the operators in on the side to the electrodeRhsFe
electrodeRhsFe->getOpPtrVector().push_back(op_loop_side);
electrodeRhsFe->getOpPtrVector().push_back(new OpElectrodeCharge<SPACE_DIM>(
CHKERR VecZeroEntries(petscVec);
"ELECTRODE", electrodeRhsFe, 0,
CHKERR VecAssemblyBegin(petscVec);
CHKERR VecAssemblyEnd(petscVec);
const double *array;
CHKERR(VecGetArrayRead(petscVec, &array));
double aLpha = array[0]; // Use explicit index instead of ZERO
double bEta = array[1]; // Use explicit index instead of ONE
MOFEM_LOG_C("SELF", Sev::inform,
"CHARGE_ELEC_1: %6.15f , CHARGE_ELEC_2: %6.15f", aLpha, bEta);
CHKERR(VecRestoreArrayRead(petscVec, &array));
}
double cal_charge_elec1;
double cal_charge_elec2;
double cal_total_energy;
const double *c_ptr, *te_ptr;
// Get a pointer to the PETSc vector data
CHKERR(VecGetArrayRead(petscVec, &c_ptr));
CHKERR(VecGetArrayRead(petscVecEnergy, &te_ptr));
// Expected charges at the electrodes
double ref_charge_elec1;
double ref_charge_elec2;
// Expected total energy of the system
double ref_tot_energy;
double tol;
cal_charge_elec1 = c_ptr[0]; // Read charge at the first electrode
cal_charge_elec2 = c_ptr[1]; // Read charge at the second electrode
cal_total_energy = te_ptr[0]; // Read total energy of the system
if (std::isnan(cal_charge_elec1) || std::isnan(cal_charge_elec2) ||
std::isnan(cal_total_energy)) {
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Atom test failed! NaN detected in calculated values.");
}
switch (atomTest) {
case 1: // 2D & 3D test
// Expected charges at the electrodes
ref_charge_elec1 = 50.0;
ref_charge_elec2 = -50.0;
// Expected total energy of the system
ref_tot_energy = 500.0;
tol = 1e-10;
break;
case 2: // wavy 3D test
ref_charge_elec1 = 10.00968352472943;
ref_charge_elec2 = 0.0; // no electrode
ref_tot_energy = 50.5978;
tol = 1e-4;
break;
default:
SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"atom test %d does not exist", atomTest);
}
// Validate the results
if (std::abs(ref_charge_elec1 - cal_charge_elec1) > tol ||
std::abs(ref_charge_elec2 - cal_charge_elec2) > tol ||
std::abs(ref_tot_energy - cal_total_energy) > tol) {
SETERRQ(
PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test %d failed! Calculated values do not match expected values",
}
CHKERR(VecRestoreArrayRead(petscVec,
&c_ptr)); // Restore the PETSc vector array
CHKERR(VecRestoreArrayRead(petscVecEnergy, &te_ptr));
}
}
//! [Get Charges]
//! [Run program]
}
//! [Run program]
//! [Main]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
Electrostatics Electrostatics_problem(m_field);
CHKERR Electrostatics_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
//! [Main]
#define MOFEM_LOG_C(channel, severity, format,...)
static char help[]
int main()
constexpr int SPACE_DIM
#define CATCH_ERRORS
Catch errors.
@ MF_ZERO
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ NOBASE
Definition definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ H1
continuous field
Definition definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ BLOCKSET
@ MOFEM_NOT_FOUND
Definition definitions.h:33
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< IntEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpInterfaceRhsVectorF
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< BASE_DIM, FIELD_DIM, SPACE_DIM > OpDomainLhsMatrixK
constexpr auto domainField
intPostProc< SPACE_DIM >::intEle IntElementForcesAndSourcesCore
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpBodySourceVectorb
const double bodySource
@ F
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition DMMoFEM.cpp:627
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:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1234
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:557
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
@ GAUSS
Gaussian quadrature integration.
@ PETSC
Standard PETSc assembly.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
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.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto bit
set bit
double D
double tol
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
auto type_from_handle(const EntityHandle h)
get type from entity handle
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
boost::shared_ptr< ForcesAndSourcesCore > electrodeRhsFe
boost::shared_ptr< ForcesAndSourcesCore > interFaceRhsFe
boost::shared_ptr< std::map< int, BlockData > > intBlockSetsPtr
SmartPetscObj< Vec > petscVec
boost::shared_ptr< DataAtIntegrationPts > commonDataPtr
MoFEM::Interface & mField
std::string domainField
MoFEMErrorCode runProgram()
[Get Charges]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
[Solve system]
MoFEMErrorCode setIntegrationRules()
[Boundary condition]
Electrostatics(MoFEM::Interface &m_field)
MoFEMErrorCode getElectrodeCharge()
[Get Total Energy]
MoFEMErrorCode readMesh()
[Read mesh]
Simple * simpleInterface
SmartPetscObj< Vec > petscVecEnergy
boost::shared_ptr< std::map< int, BlockData > > permBlockSetsPtr
MoFEMErrorCode assembleSystem()
[Set integration rules]
PetscBool is_partitioned
boost::shared_ptr< std::map< int, BlockData > > electrodeBlockSetsPtr
MoFEMErrorCode getTotalEnergy()
[Output results]
MoFEMErrorCode boundaryCondition()
[Setup problem]
MoFEMErrorCode solveSystem()
[Assemble system]
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
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:118
Deprecated interface functions.
MoFEMErrorCode buildProblem()
Build problem.
Definition Simple.cpp:722
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:261
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition Simple.cpp:471
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition Simple.cpp:660
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:355
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800
MoFEMErrorCode buildFields()
Build fields.
Definition Simple.cpp:584
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition Simple.cpp:551
MoFEMErrorCode addDataField(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 data field.
Definition Simple.cpp:393
bool & getAddSkeletonFE()
Get the addSkeletonFE flag.
Definition Simple.hpp:536
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:736
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:450
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:429
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.