28 using namespace MoFEM;
30 static char help[] =
"...\n\n";
68 :
public boost::enable_shared_from_this<Example::CommonData> {
73 return boost::shared_ptr<VectorDouble>(shared_from_this(),
74 &rhoAtIntegrationPts);
102 OpZero(boost::shared_ptr<CommonData> &common_data_ptr)
103 :
OpElement(
"rho", OPROW), commonDataPtr(common_data_ptr) {
104 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
113 OpFirst(boost::shared_ptr<CommonData> &common_data_ptr)
114 :
OpElement(
"rho", OPROW), commonDataPtr(common_data_ptr) {
115 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
125 OpSecond(boost::shared_ptr<CommonData> &common_data_ptr)
126 :
OpElement(
"rho", OPROW), commonDataPtr(common_data_ptr) {
127 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
140 CHKERR createCommonData();
143 CHKERR integrateElements();
158 constexpr
int order = 1;
168 commonDataPtr = boost::make_shared<CommonData>();
171 if (mField.get_comm_rank() == 0)
173 local_size = CommonData::LAST_ELEMENT;
179 CommonData::LAST_ELEMENT);
188 auto set_density = [&](
VectorAdaptor &&field_data,
double *xcoord,
189 double *ycoord,
double *zcoord) {
195 CHKERR mField.getInterface(field_blas);
209 "rho", commonDataPtr->getRhoAtIntegrationPtsPtr()));
235 CHKERR VecZeroEntries(commonDataPtr->petscVec);
242 CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
243 CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
252 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
253 if (mField.get_comm_rank() == 0) {
254 MOFEM_LOG_C(
"SELF", Sev::inform,
"Mass %6.4e", array[CommonData::ZERO]);
256 "First moment of inertia [ %6.4e, %6.4e, %6.4e ]",
257 array[CommonData::FIRST_X], array[CommonData::FIRST_Y],
258 array[CommonData::FIRST_Z]);
260 "Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
262 array[CommonData::SECOND_XX], array[CommonData::SECOND_XY],
263 array[CommonData::SECOND_XZ], array[CommonData::SECOND_YY],
264 array[CommonData::SECOND_YZ], array[CommonData::SECOND_ZZ]);
266 CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
275 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
277 PetscBool test = PETSC_FALSE;
279 if (mField.get_comm_rank() == 0 && test) {
280 constexpr
double eps = 1e-8;
281 constexpr
double expected_volume = 1.;
282 constexpr
double expected_first_moment = 0.;
283 constexpr
double expected_second_moment = 1. / 12.;
284 if (std::abs(array[CommonData::ZERO] - expected_volume) >
eps)
286 "Wrong area %6.4e != %6.4e", expected_volume,
287 array[CommonData::ZERO]);
289 {CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z}) {
290 if (std::abs(array[
i] - expected_first_moment) >
eps)
292 "Wrong first moment %6.4e != %6.4e", expected_first_moment,
295 for (
auto i : {CommonData::SECOND_XX, CommonData::SECOND_YY,
296 CommonData::SECOND_ZZ}) {
297 if (std::abs(array[
i] - expected_second_moment) >
eps)
299 "Wrong second moment %6.4e != %6.4e", expected_second_moment,
303 CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
310 int main(
int argc,
char *argv[]) {
320 DMType dm_name =
"DMMOFEM";
347 if (
type == MBVERTEX) {
349 const int nb_integration_pts =
350 getGaussPts().size2();
351 auto t_w = getFTensor0IntegrationWeight();
353 commonDataPtr->rhoAtIntegrationPts);
354 const double volume = getMeasure();
357 double element_local_value = 0;
358 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
359 element_local_value += t_w * t_rho * volume;
365 const int index = CommonData::ZERO;
366 CHKERR VecSetValue(commonDataPtr->petscVec,
index, element_local_value,
377 const int nb_integration_pts = getGaussPts().size2();
378 auto t_w = getFTensor0IntegrationWeight();
380 commonDataPtr->rhoAtIntegrationPts);
382 getFTensor1CoordsAtGaussPts();
383 const double volume = getMeasure();
392 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
394 t_s(
i) += t_w * t_rho * volume * t_x(
i);
402 constexpr std::array<int, 3> indices = {
403 CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z};
417 const int nb_integration_pts = getGaussPts().size2();
418 auto t_w = getFTensor0IntegrationWeight();
420 auto t_x = getFTensor1CoordsAtGaussPts();
421 const double volume = getMeasure();
424 std::array<double, 6> element_local_value;
425 std::fill(element_local_value.begin(), element_local_value.end(), 0.0);
429 &element_local_value[0], &element_local_value[1], &element_local_value[2],
430 &element_local_value[3], &element_local_value[4],
431 &element_local_value[5]);
437 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
440 t_I(
i,
j) += (t_w * t_rho * volume) * (t_x(
i) ^ t_x(
j));
448 constexpr std::array<int, 6> indices = {
449 CommonData::SECOND_XX, CommonData::SECOND_XY, CommonData::SECOND_XZ,
450 CommonData::SECOND_YY, CommonData::SECOND_YZ, CommonData::SECOND_ZZ};
454 &element_local_value[0], ADD_VALUES);
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
#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
#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.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
FTensor::Index< 'i', SPACE_DIM > i
int main(int argc, char *argv[])
[Test example]
EntitiesFieldData::EntData EntData
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< double > VectorAdaptor
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
auto createSmartVectorMPI
Create MPI Vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::shared_ptr< VectorDouble > getRhoAtIntegrationPtsPtr()
VecElements
Vector to indicate indices for storing zero, first and second moments of inertia.
SmartPetscObj< Vec > petscVec
Smart pointer which stores PETSc distributed vector.
VectorDouble rhoAtIntegrationPts
Storing density at integration points.
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[ZeroOp]
OpFirst(boost::shared_ptr< CommonData > &common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[FirstOp]
boost::shared_ptr< CommonData > commonDataPtr
OpSecond(boost::shared_ptr< CommonData > &common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[main]
boost::shared_ptr< CommonData > commonDataPtr
OpZero(boost::shared_ptr< CommonData > &common_data_ptr)
MoFEMErrorCode pushOperators()
[Set density distribution]
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode setUp()
[Run all]
MoFEMErrorCode checkResults()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEMErrorCode postProcess()
MoFEMErrorCode setFieldValues()
[Create common data]
MoFEMErrorCode integrateElements()
[Push operators to pipeline]
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)
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
Get value at integration points for scalar field.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
friend class UserDataOperator