16static char help[] =
"...\n\n";
18struct OpVolume :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
29 const int nb_int_pts = getGaussPts().size2();
31 auto t_w = getFTensor0IntegrationWeight();
32 double v = getMeasure();
34 for (
int gg = 0; gg != nb_int_pts; gg++) {
38 CHKERR VecSetValue(
vOl, 0, vol, ADD_VALUES);
51struct OpFace :
public FaceElementForcesAndSourcesCore::UserDataOperator {
68 for (
int gg = 0; gg != nb_int_pts; gg++) {
69 vol += (t_coords(
i) * t_normal(
i)) * t_w;
75 CHKERR VecSetValue(
vOl, 0, vol, ADD_VALUES);
95int main(
int argc,
char *argv[]) {
103 moab::Core moab_core;
104 moab::Interface &moab = moab_core;
111 DMType dm_name =
"DMMOFEM";
138 auto dm = simple_interface->
getDM();
141 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
143 boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
146 domain_fe->getRuleHook =
VolRule();
147 boundary_fe->getRuleHook =
FaceRule();
157 auto material_grad_mat = boost::make_shared<MatrixDouble>();
158 auto material_det_vec = boost::make_shared<VectorDouble>();
160 domain_fe->getOpPtrVector().push_back(
164 material_grad_mat, material_det_vec,
nullptr));
165 domain_fe->getOpPtrVector().push_back(
167 domain_fe->getOpPtrVector().push_back(
169 domain_fe->getOpPtrVector().push_back(
170 new OpVolume(
"MESH_NODE_POSITIONS", vol));
172 boundary_fe->getOpPtrVector().push_back(
174 boundary_fe->getOpPtrVector().push_back(
175 new OpFace(
"MESH_NODE_POSITIONS", surf_vol));
183 auto skeleton_fe = boost::make_shared<FEMethod>();
194 skeleton_fe->x_t = x_t;
195 skeleton_fe->x_tt = x_tt;
197 skeleton_fe->preProcessHook = [&]() {
199 if (f != skeleton_fe->ts_F)
201 if (A != skeleton_fe->ts_A)
203 if (
B != skeleton_fe->ts_B)
205 if (x != skeleton_fe->ts_u)
207 if (x_t != skeleton_fe->ts_u_t)
209 if (x_tt != skeleton_fe->ts_u_tt)
214 skeleton_fe->postProcessHook = []() {
return 0; };
215 skeleton_fe->operatorHook = []() {
return 0; };
221 CHKERR VecAssemblyBegin(vol);
222 CHKERR VecAssemblyEnd(vol);
223 CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
224 CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
225 CHKERR VecAssemblyBegin(surf_vol);
226 CHKERR VecAssemblyEnd(surf_vol);
227 CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
228 CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
231 CHKERR VecGetArray(vol, &a_vol);
233 CHKERR VecGetArray(surf_vol, &a_surf_vol);
234 cout <<
"Volume = " << a_vol[0] << endl;
235 cout <<
"Surf Volume = " << a_surf_vol[0] << endl;
236 if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
239 CHKERR VecRestoreArray(vol, &a_vol);
240 CHKERR VecRestoreArray(vol, &a_surf_vol);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#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.
#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.
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.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
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
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 auto field_name
Set integration rule to boundary elements.
int operator()(int, int, int) const
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor1NormalsAtGaussPts()
get normal at integration points
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.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Set inverse jacobian to base functions.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
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 getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
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 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.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
OpFace(const std::string &field_name, Vec vol)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpVolume(const std::string &field_name, Vec vol)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
int operator()(int, int, int) const