#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
static char help[] =
"...\n\n";
public:
private:
};
true;
}
true);
auto get_ents_by_dim = [&](const auto dim) {
return domain_ents;
} else {
if (dim == 0)
else
return ents;
}
};
auto get_base = [&]() {
auto domain_ents = get_ents_by_dim(
SPACE_DIM);
if (domain_ents.empty())
switch (type) {
case MBQUAD:
case MBHEX:
case MBTRI:
case MBTET:
default:
}
};
auto base = get_base();
PETSC_NULLPTR);
PETSC_NULLPTR);
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(
mField,
"GEOMETRY");
};
if (
bit->getName().compare(0, 12,
"MAT_ELECTRIC") == 0) {
const int id =
bit->getMeshsetId();
auto &block_data = (*permBlockSetsPtr)[id];
mat_electr_ents.merge(block_data.domainEnts);
std::vector<double> attributes;
bit->getAttributes(attributes);
if (attributes.size() < 1) {
" At least one permittivity attributes should be given but "
"found %d",
attributes.size());
}
block_data.iD = id;
block_data.epsPermit = attributes[0];
}
}
if (
bit->getName().compare(0, 12,
"INT_ELECTRIC") == 0) {
const int id =
bit->getMeshsetId();
auto &block_data = (*intBlockSetsPtr)[id];
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) {
"At least one charge attributes should be given but found %d",
attributes.size());
}
block_data.iD = id;
block_data.chargeDensity = attributes[0];
}
}
int electrodeCount = 0;
if (
bit->getName().compare(0, 9,
"ELECTRODE") == 0) {
const int id =
bit->getMeshsetId();
auto &block_data = (*electrodeBlockSetsPtr)[id];
++electrodeCount;
bit->getMeshset(),
SPACE_DIM - 1, block_data.electrodeEnts,
true);
electrode_ents.merge(block_data.electrodeEnts);
if (electrodeCount > 2) {
"Three or more electrode blocksets found");
;
}
}
}
mat_electr_ents);
int_electr_ents);
electrode_ents);
PETSC_NULLPTR);
PETSC_NULLPTR);
CHKERR skinner.find_skin(0, mat_electr_ents,
false, skin_tris);
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin_tris,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &proc_skin);
} else {
proc_skin = skin_tris;
}
"SKIN");
"ELECTRODE");
DMType dm_name = "DMMOFEM";
SmartPetscObj<DM> dm;
CHKERR DMSetType(dm, dm_name);
int local_size;
else
local_size = 0;
}
CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
}
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; };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
}
auto add_domain_base_ops = [&](auto &pipeline) {
"GEOMETRY");
pipeline.push_back(
};
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
};
{
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{
auto set_values_to_bc_dofs = [&](auto &fe) {
auto get_bc_hook = [&]() {
EssentialPreProc<TemperatureCubitBcData> hook(
mField, fe, {});
return hook;
};
fe->preProcessHook = get_bc_hook();
};
auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
using OpInternal =
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));
};
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());
};
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
};
{
};
}
}
auto ksp_solver = pipeline_mng->createKSP();
boost::shared_ptr<ForcesAndSourcesCore> null;
DM dm;
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
double fnorm;
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"F norm = %9.8e\n", fnorm);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
double dnorm;
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"D norm = %9.8e\n", dnorm);
}
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);
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>();
"GEOMETRY");
pipeline.push_back(
new OpCalculateScalarFieldValues(
domainField, u_ptr));
pipeline.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
pipeline.push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>(
domainField, 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(
post_proc_fe->getOpPtrVector().push_back(
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");
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(
permBlockSetsPtr, commonDataPtr));
op_loop_skin->getOpPtrVector().push_back(
permBlockSetsPtr, commonDataPtr));
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{}));
post_proc_skin);
CHKERR post_proc_skin->writeFile(
"out_skin.h5m");
}
}
pip_energy->getDomainRhsFE().reset();
pip_energy->getDomainLhsFE().reset();
pip_energy->getBoundaryLhsFE().reset();
pip_energy->getBoundaryRhsFE().reset();
pip_energy->getOpDomainRhsPipeline().clear();
pip_energy->getOpDomainLhsPipeline().clear();
boost::shared_ptr<std::map<int, BlockData>> intrnlDomnBlckSetPtr =
boost::make_shared<std::map<int, BlockData>>();
if (
bit->getName().compare(0, 10,
"DOMAIN_INT") == 0) {
const int id =
bit->getMeshsetId();
auto &block_data = (*intrnlDomnBlckSetPtr)[id];
bit->getMeshset(),
SPACE_DIM, block_data.internalDomainEnts,
true);
internal_domain.merge(block_data.internalDomainEnts);
block_data.iD = id;
}
}
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(
pip_energy->getOpDomainRhsPipeline().push_back(
pip_energy->loopFiniteElements();
double total_energy = 0.0;
const double *array;
total_energy = array[
ZERO];
MOFEM_LOG_C(
"SELF", Sev::inform,
"Total Energy: %6.15f", total_energy);
}
}
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(
auto d_jump = boost::make_shared<MatrixDouble>();
};
const double *array;
"CHARGE_ELEC_1: %6.15f , CHARGE_ELEC_2: %6.15f",
aLpha,
bEta);
}
double cal_charge_elec1;
double cal_charge_elec2;
double cal_total_energy;
const double *c_ptr, *te_ptr;
double ref_charge_elec1;
double ref_charge_elec2;
double ref_tot_energy;
cal_charge_elec1 = c_ptr[0];
cal_charge_elec2 = c_ptr[1];
cal_total_energy = te_ptr[0];
if (std::isnan(cal_charge_elec1) || std::isnan(cal_charge_elec2) ||
std::isnan(cal_total_energy)) {
"Atom test failed! NaN detected in calculated values.");
}
case 1:
ref_charge_elec1 = 50.0;
ref_charge_elec2 = -50.0;
ref_tot_energy = 500.0;
break;
case 2:
ref_charge_elec1 = 10.00968352472943;
ref_charge_elec2 = 0.0;
ref_tot_energy = 50.5978;
break;
default:
"atom test %d does not exist",
atomTest);
}
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(
"atom test %d failed! Calculated values do not match expected values",
}
&c_ptr));
}
}
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
CHKERR Electrostatics_problem.runProgram();
}
return 0;
}
#define MOFEM_LOG_C(channel, severity, format,...)
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#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 ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#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
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
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
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
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.
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
#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.
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
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]
SmartPetscObj< Vec > petscVecEnergy
boost::shared_ptr< std::map< int, BlockData > > permBlockSetsPtr
MoFEMErrorCode assembleSystem()
[Set integration rules]
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
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
MoFEMErrorCode buildProblem()
Build problem.
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.
MoFEMErrorCode defineFiniteElements()
Define finite elements.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode buildFiniteElements()
Build finite elements.
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.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode buildFields()
Build fields.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
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.
bool & getAddSkeletonFE()
Get the addSkeletonFE flag.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getProblemName() const
Get the Problem Name.
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.