25 using namespace MoFEM;
27 static char help[] =
"...\n\n";
53 static double fUn(
const double x,
const double y,
double z) {
56 for (
int i = 0;
i <=
o; ++
i) {
59 r += pow(x,
i) * pow(y,
j);
69 for (
int i = 0;
i <=
o; ++
i) {
72 r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) : 0;
73 r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) : 0;
82 static double fUn(
const double x,
const double y,
double z) {
85 for (
int i = 0;
i <=
o; ++
i) {
86 for (
int j = 0;
j <=
o -
i;
j++) {
89 r += pow(x,
i) * pow(y,
j) * pow(z,
k);
101 for (
int i = 0;
i <=
o; ++
i) {
102 for (
int j = 0;
j <=
o -
i;
j++) {
105 r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) * pow(z,
k) : 0;
106 r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) * pow(z,
k) : 0;
107 r(2) +=
k > 0 ?
k * pow(x,
i) * pow(y,
j) * pow(z,
k - 1) : 0;
123 :
DomainEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
124 checkGradients(check_grads) {}
131 const int nb_gauss_pts = getGaussPts().size2();
132 if (
type == MBVERTEX) {
133 vAls.resize(nb_gauss_pts,
false);
134 diffVals.resize(
SPACE_DIM, nb_gauss_pts,
false);
143 <<
"Type " << moab::CN::EntityTypeName(
type) <<
" side " << side;
149 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
151 for (
int bb = 0; bb != nb_dofs; bb++) {
152 t_vals += t_base_fun * t_data;
156 const double v = t_vals;
157 if (!std::isnormal(
v))
162 if (checkGradients) {
163 auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
165 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
167 for (
int bb = 0; bb != nb_dofs; bb++) {
168 t_diff_vals(
i) += t_diff_base_fun(
i) * t_data;
173 if (!std::isnormal(t_diff_vals(
d)))
191 boost::shared_ptr<VectorDouble> &ptr_vals,
192 boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
194 :
DomainEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
195 ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals),
196 checkGradients(check_grads) {
197 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
205 const double eps = 1e-6;
206 const int nb_gauss_pts = data.
getN().size1();
212 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
217 err_val = std::abs(t_vals - t_ptr_vals);
218 MOFEM_LOG(
"AT", Sev::noisy) <<
"Val op error " << err_val;
222 "Wrong value from operator %4.3e", err_val);
224 const double x = getCoordsAtGaussPts()(gg, 0);
225 const double y = getCoordsAtGaussPts()(gg, 1);
226 const double z = getCoordsAtGaussPts()(gg, 2);
231 err_val = std::fabs(delta_val * delta_val);
232 MOFEM_LOG(
"AT", Sev::verbose) << err_val <<
" : " << t_vals;
241 if (checkGradients) {
243 auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
244 auto t_ptr_diff_vals = getFTensor1FromMat<SPACE_DIM>(*ptrDiffVals);
246 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
251 t_delta_diff_val(
i) = t_diff_vals(
i) - t_ptr_diff_vals(
i);
252 err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
253 MOFEM_LOG(
"AT", Sev::noisy) <<
"Diff val op error " << err_diff_val;
255 if (err_diff_val >
eps)
257 "Wrong derivatives from operator %4.3e", err_diff_val);
259 const double x = getCoordsAtGaussPts()(gg, 0);
260 const double y = getCoordsAtGaussPts()(gg, 1);
261 const double z = getCoordsAtGaussPts()(gg, 2);
265 t_delta_diff_val(
i) = t_diff_vals(
i) - t_diff_anal(
i);
267 err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
270 <<
"Diff val " << err_diff_val <<
" : "
271 << sqrt(t_diff_vals(
i) * t_diff_vals(
i)) <<
" : "
272 << t_diff_vals(0) <<
" (" << t_diff_anal(0) <<
") "
273 << t_diff_vals(1) <<
" (" << t_diff_anal(1) <<
") "
274 << t_diff_vals(2) <<
" (" << t_diff_anal(2) <<
")";
277 <<
"Diff val " << err_diff_val <<
" : "
278 << sqrt(t_diff_vals(
i) * t_diff_vals(
i)) <<
" : "
279 << t_diff_vals(0) <<
" (" << t_diff_anal(0) <<
") "
280 << t_diff_vals(1) <<
" (" << t_diff_anal(1) <<
")";
283 << getCoords()(3 * 1 + 0) - getCoords()(3 * 0 + 0);
285 << getCoords()(3 * 1 + 1) - getCoords()(3 * 0 + 1);
287 << getCoords()(3 * 1 + 2) - getCoords()(3 * 0 + 2);
289 MOFEM_LOG(
"AT", Sev::verbose) <<
"Diff val error " << err_diff_val;
290 if (err_diff_val >
eps)
292 "Wrong derivative of value %4.3e %4.3e", err_diff_val,
304 int main(
int argc,
char *argv[]) {
310 DMType dm_name =
"DMMOFEM";
317 auto core_log = logging::core::get();
340 const char *list_bases[] = {
"ainsworth",
"ainsworth_labatto",
"demkowicz",
343 PetscInt choice_base_value = AINSWORTH;
345 LASBASETOP, &choice_base_value, &flg);
347 if (flg != PETSC_TRUE)
350 if (choice_base_value == AINSWORTH)
352 if (choice_base_value == AINSWORTH_LOBATTO)
354 else if (choice_base_value == DEMKOWICZ)
356 else if (choice_base_value == BERNSTEIN)
359 enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
360 const char *list_spaces[] = {
"h1",
"l2"};
361 PetscInt choice_space_value = H1SPACE;
363 LASBASETSPACE, &choice_space_value, &flg);
364 if (flg != PETSC_TRUE)
367 if (choice_space_value == H1SPACE)
369 else if (choice_space_value == L2SPACE)
379 CHKERR moab.get_entities_by_dimension(0, 1, edges);
380 CHKERR moab.get_entities_by_dimension(0, 2, faces);
382 if (choice_base_value != BERNSTEIN) {
386 for (
auto e : edges) {
388 rise_order.insert(e);
393 for (
auto f : faces) {
395 rise_order.insert(
f);
401 for (
auto e : rise_order) {
403 rise_twice.insert(e);
416 auto dm = simple_interface->
getDM();
419 auto jac_ptr = boost::make_shared<MatrixDouble>();
420 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
421 auto det_ptr = boost::make_shared<VectorDouble>();
424 auto assemble_matrices_and_vectors = [&]() {
427 pipeline_mng->getOpDomainRhsPipeline().push_back(
434 pipeline_mng->getOpDomainLhsPipeline().push_back(
436 pipeline_mng->getOpDomainLhsPipeline().push_back(
438 pipeline_mng->getOpDomainLhsPipeline().push_back(
443 pipeline_mng->getOpDomainLhsPipeline().push_back(
445 pipeline_mng->getOpDomainLhsPipeline().push_back(
447 pipeline_mng->getOpDomainLhsPipeline().push_back(
449 pipeline_mng->getOpDomainRhsPipeline().push_back(
451 pipeline_mng->getOpDomainRhsPipeline().push_back(
453 pipeline_mng->getOpDomainRhsPipeline().push_back(
459 pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass(
460 "FIELD1",
"FIELD1", [](
double,
double,
double) {
return 1.; }));
461 pipeline_mng->getOpDomainRhsPipeline().push_back(
465 return 2 * p_data + 1;
473 auto solve_problem = [&] {
475 auto solver = pipeline_mng->createKSP();
476 CHKERR KSPSetFromOptions(solver);
479 auto dm = simple_interface->
getDM();
484 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
485 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
491 auto check_solution = [&] {
494 boost::shared_ptr<VectorDouble> ptr_values(
new VectorDouble());
495 boost::shared_ptr<MatrixDouble> ptr_diff_vals(
new MatrixDouble());
497 pipeline_mng->getDomainLhsFE().reset();
498 pipeline_mng->getOpDomainRhsPipeline().clear();
501 pipeline_mng->getOpDomainRhsPipeline().push_back(
503 pipeline_mng->getOpDomainRhsPipeline().push_back(
505 pipeline_mng->getOpDomainRhsPipeline().push_back(
510 pipeline_mng->getOpDomainRhsPipeline().push_back(
512 pipeline_mng->getOpDomainRhsPipeline().push_back(
514 pipeline_mng->getOpDomainRhsPipeline().push_back(
518 pipeline_mng->getOpDomainRhsPipeline().push_back(
520 pipeline_mng->getOpDomainRhsPipeline().push_back(
522 pipeline_mng->getOpDomainRhsPipeline().push_back(
526 vals, diff_vals, ptr_values, ptr_diff_vals,
true));
528 CHKERR pipeline_mng->loopFiniteElements();
533 CHKERR assemble_matrices_and_vectors();
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
#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.
auto smartCreateDMVector
Get smart vector from 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.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
const double r
rate factor
#define EXECUTABLE_DIMENSION
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
static FTensor::Tensor1< double, 3 > fUn(const double x, const double y)
static FTensor::Tensor2< double, 3, 2 > diffFun(const double x, const double y)
static double fUn(const double x, const double y, double z)
static FTensor::Tensor1< double, 2 > diffFun(const double x, const double y, double z)
static FTensor::Tensor1< double, 3 > diffFun(const double x, const double y, double z)
static double fUn(const double x, const double y, double z)
PipelineManager::FaceEle DomainEle
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)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Calculate jacobian on Hex or other volume which is not simplex.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
Transform local reference derivatives of shape functions to global derivatives.
PipelineManager interface.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Simple interface for fast problem set-up.
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 setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element with switches.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
OpCheckValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, boost::shared_ptr< VectorDouble > &ptr_vals, boost::shared_ptr< MatrixDouble > &ptr_diff_vals, bool check_grads)
boost::shared_ptr< VectorDouble > ptrVals
const bool checkGradients
FTensor::Index< 'i', SPACE_DIM > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > ptrDiffVals
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const bool checkGradients
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
FTensor::Index< 'i', SPACE_DIM > i