16static char help[] =
"...\n\n";
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);
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->meshPositionsFieldName =
"none";
161 domain_fe->getOpPtrVector().push_back(
165 material_grad_mat, material_det_vec,
nullptr));
166 domain_fe->getOpPtrVector().push_back(
168 domain_fe->getOpPtrVector().push_back(
170 domain_fe->getOpPtrVector().push_back(
171 new OpVolume(
"MESH_NODE_POSITIONS", vol));
173 boundary_fe->getOpPtrVector().push_back(
175 boundary_fe->getOpPtrVector().push_back(
176 new OpFace(
"MESH_NODE_POSITIONS", surf_vol));
184 auto skeleton_fe = boost::make_shared<FEMethod>();
195 skeleton_fe->x_t = x_t;
196 skeleton_fe->x_tt = x_tt;
198 skeleton_fe->preProcessHook = [&]() {
200 if (f != skeleton_fe->ts_F)
202 if (
A != skeleton_fe->ts_A)
204 if (
B != skeleton_fe->ts_B)
206 if (x != skeleton_fe->ts_u)
208 if (x_t != skeleton_fe->ts_u_t)
210 if (x_tt != skeleton_fe->ts_u_tt)
215 skeleton_fe->postProcessHook = []() {
return 0; };
216 skeleton_fe->operatorHook = []() {
return 0; };
222 CHKERR VecAssemblyBegin(vol);
223 CHKERR VecAssemblyEnd(vol);
224 CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
225 CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
226 CHKERR VecAssemblyBegin(surf_vol);
227 CHKERR VecAssemblyEnd(surf_vol);
228 CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
229 CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
232 CHKERR VecGetArray(vol, &a_vol);
234 CHKERR VecGetArray(surf_vol, &a_surf_vol);
235 cout <<
"Volume = " << a_vol[0] << endl;
236 cout <<
"Surf Volume = " << a_surf_vol[0] << endl;
237 if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
240 CHKERR VecRestoreArray(vol, &a_vol);
241 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)
default operator for TRI element
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.
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 loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
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