1#ifndef __STDOPERATORS_HPP__
2#define __STDOPERATORS_HPP__
10using OpEle = FaceElementForcesAndSourcesCore::UserDataOperator;
15using EntData = DataForcesAndSourcesCore::EntData;
42 coef.resize(3, 3,
false);
43 rate.resize(3,
false);
56 for (
int i = 0;
i < 3; ++
i) {
85 const int nb_row_dofs = row_data.getIndices().size();
86 const int nb_col_dofs = col_data.getIndices().size();
87 if (nb_row_dofs && nb_col_dofs) {
88 const int nb_integration_pts =
getGaussPts().size2();
89 mat.resize(nb_row_dofs, nb_col_dofs,
false);
91 auto t_row_base = row_data.getFTensor0N();
94 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
95 const double a = t_w * vol;
96 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
97 auto t_col_base = col_data.getFTensor0N(gg, 0);
98 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
99 mat(rr, cc) +=
a * t_row_base * t_col_base;
107 if (row_side != col_side || row_type != col_type) {
108 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
123 boost::shared_ptr<PreviousData> &data1,
124 boost::shared_ptr<PreviousData> &data2,
125 boost::shared_ptr<PreviousData> &data3,
126 std::map<int, BlockData> &block_map)
131 boost::shared_ptr<VectorDouble> slow_value_ptr1(
commonData1,
133 boost::shared_ptr<VectorDouble> slow_value_ptr2(
commonData2,
135 boost::shared_ptr<VectorDouble> slow_value_ptr3(
commonData3,
141 const int nb_integration_pts =
getGaussPts().size2();
142 if (type == MBVERTEX) {
143 vec1.resize(nb_integration_pts,
false);
144 vec2.resize(nb_integration_pts,
false);
145 vec3.resize(nb_integration_pts,
false);
150 const int nb_dofs = data.getIndices().size();
153 auto find_block_data = [&]() {
157 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
158 block_raw_ptr = &
m.second;
162 return block_raw_ptr;
165 auto block_data_ptr = find_block_data();
169 auto &block_data = *block_data_ptr;
171 const int nb_integration_pts =
getGaussPts().size2();
183 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
184 t_slow_values1 = block_data.rate[0] * t_mass_values1 *
185 (1.0 - block_data.coef(0, 0) * t_mass_values1 -
186 block_data.coef(0, 1) * t_mass_values2 -
187 block_data.coef(0, 2) * t_mass_values3);
188 t_slow_values2 = block_data.rate[1] * t_mass_values2 *
189 (1.0 - block_data.coef(1, 0) * t_mass_values1 -
190 block_data.coef(1, 1) * t_mass_values2 -
191 block_data.coef(1, 2) * t_mass_values3);
193 t_slow_values3 = block_data.rate[2] * t_mass_values3 *
194 (1.0 - block_data.coef(2, 0) * t_mass_values1 -
195 block_data.coef(2, 1) * t_mass_values2 -
196 block_data.coef(2, 2) * t_mass_values3);
218 typedef boost::function<
double(
const double,
const double,
const double)>
220 typedef boost::function<FTensor::Tensor1<double, 3>(
221 const double,
const double,
const double)>
224 boost::shared_ptr<PreviousData> &data,
234 const int nb_dofs = data.getIndices().size();
236 vecF.resize(nb_dofs,
false);
238 const int nb_integration_pts =
getGaussPts().size2();
242 auto t_base = data.getFTensor0N();
249 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
250 const double a = vol * t_w;
252 double u_dot =
exactDot(t_coords(
NX), t_coords(
NY), ct);
253 double u_lap =
exactLap(t_coords(
NX), t_coords(
NY), ct);
256 const double f = t_slow_value;
258 for (
int rr = 0; rr != nb_dofs; ++rr) {
259 const double b =
a * f * t_base;
290 typedef boost::function<
double(
const double,
const double,
const double)>
292 typedef boost::function<FTensor::Tensor1<double, 3>(
293 const double,
const double,
const double)>
297 std::map<int, BlockData> &block_map,
FVal exact_value,
304 const int nb_dofs = data.getIndices().size();
306 auto find_block_data = [&]() {
310 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
311 block_raw_ptr = &
m.second;
315 return block_raw_ptr;
318 auto block_data_ptr = find_block_data();
322 auto &block_data = *block_data_ptr;
324 vecF.resize(nb_dofs,
false);
326 const int nb_integration_pts =
getGaussPts().size2();
330 auto t_grad = getFTensor1FromMat<DIM>(
commonData->grads);
332 auto t_base = data.getFTensor0N();
333 auto t_diff_base = data.getFTensor1DiffN<DIM>();
345 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
346 const double a = vol * t_w;
348 double u_dot =
exactDot(t_coords(
NX), t_coords(
NY), ct);
353 -block_data.B0 *
exactLap(t_coords(
NX), t_coords(
NY), ct);
358 for (
int rr = 0; rr != nb_dofs; ++rr) {
360 a * (t_base * t_dot_val +
361 (block_data.B0 +
B * t_val) * t_diff_base(
i) * t_grad(
i) -
398 std::map<int, BlockData> &block_map)
408 const int nb_row_dofs = row_data.getIndices().size();
409 const int nb_col_dofs = col_data.getIndices().size();
412 if (nb_row_dofs && nb_col_dofs) {
413 auto find_block_data = [&]() {
417 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
418 block_raw_ptr = &
m.second;
422 return block_raw_ptr;
425 auto block_data_ptr = find_block_data();
429 auto &block_data = *block_data_ptr;
431 mat.resize(nb_row_dofs, nb_col_dofs,
false);
433 const int nb_integration_pts =
getGaussPts().size2();
434 auto t_row_base = row_data.getFTensor0N();
437 auto t_grad = getFTensor1FromMat<DIM>(
commonData->grads);
439 auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
448 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
449 const double a = vol * t_w;
451 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
452 auto t_col_base = col_data.getFTensor0N(gg, 0);
453 auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
454 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
457 a * (t_row_base * t_col_base * ts_a +
458 (block_data.B0 +
B * t_val) * t_row_diff_base(
i) *
460 B * t_col_base * t_grad(
i) * t_row_diff_base(
i));
474 if (row_side != col_side || row_type != col_type) {
475 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
495 cerr <<
"OpAssembleNaturalBCRhsTau()" << endl;
500 const int nb_dofs = data.getIndices().size();
508 vecF.resize(nb_dofs,
false);
510 const int nb_integration_pts =
getGaussPts().size2();
511 auto t_row_base = data.getFTensor0N();
517 const double pi = 3.141592653589793;
519 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
520 const double a = t_w;
524 double mm = -10 * 8 * pi * cos(2 * pi * x) * sin(2 * pi * y) *
526 double nn = -10 * 8 * pi * sin(2 * pi * x) * cos(2 * pi * y) *
530 double h = t_bd_val(
i) * t_normal(
i);
531 for (
int rr = 0; rr != nb_dofs; ++rr) {
532 vecF[rr] += t_row_base *
h *
a;
552 typedef boost::function<
double(
const double,
const double,
const double)>
554 typedef boost::function<FTensor::Tensor1<double, 3>(
555 const double,
const double,
const double)>
559 boost::shared_ptr<PreviousData> &prev_data,
560 std::map<int, BlockData> &block_map,
double &err)
566 const int nb_dofs = data.getFieldData().size();
570 auto find_block_data = [&]() {
574 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
575 block_raw_ptr = &
m.second;
579 return block_raw_ptr;
582 auto block_data_ptr = find_block_data();
586 auto &block_data = *block_data_ptr;
589 auto t_grad = getFTensor1FromMat<2>(
prevData->grads);
592 data.getFieldData().clear();
594 const int nb_integration_pts =
getGaussPts().size2();
604 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
605 const double a = vol * t_w;
607 double mass_exact =
exactVal(t_coords(
NX), t_coords(
NY), ct);
609 -block_data.B0 *
exactLap(t_coords(
NX), t_coords(
NY), ct);
617 pow(block_data.B0, 2) * (pow(flux_exact(0) - t_grad(0), 2) +
618 pow(flux_exact(1) - t_grad(1), 2));
620 double local_error = pow(mass_exact - t_value, 2);
622 data.getFieldData()[0] +=
a * local_error;
631 data.getFieldDofs()[0]->getFieldData() = data.getFieldData()[0];
650 boost::shared_ptr<PostProcFaceOnRefinedMesh> &post_proc,
double &err)
661 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
664 CHKERR VecCreateMPI(
cOmm, 1, PETSC_DECIDE, &error_per_proc);
665 auto get_global_error = [&]() {
670 CHKERR get_global_error();
671 CHKERR VecAssemblyBegin(error_per_proc);
672 CHKERR VecAssemblyEnd(error_per_proc);
674 CHKERR VecSum(error_per_proc, &error_sum);
675 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Error : %3.4e \n", error_sum);
682 boost::shared_ptr<PostProcFaceOnRefinedMesh>
postProc;
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
FTensor::Index< 'i', 3 > i
DataForcesAndSourcesCore::EntData EntData
const double essential_bc_values
const double natural_bc_values
bool sYmm
If true assume that matrix is symmetric structure.
default operator for EDGE element
VectorDouble & getDirection()
get edge direction
structure for User Loop Methods on finite elements
default operator for TRI element
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
intrusive_ptr for managing petsc objects
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProc
MoFEMErrorCode operator()()
function is run for every finite element
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcFaceOnRefinedMesh > &post_proc, double &err)
MoFEMErrorCode postProcess()
function is run at the end of loop
OpAssembleMass(std::string fieldu, SmartPetscObj< Mat > m)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleNaturalBCRhs(std::string mass_field, Range &natural_bd_ents)
boost::shared_ptr< PreviousData > commonData
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleSlowRhs(std::string field, boost::shared_ptr< PreviousData > &data, FVal exact_value, FVal exact_dot, FVal exact_lap)
boost::function< double(const double, const double, const double)> FVal
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleStiffLhs(std::string fieldu, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
OpAssembleStiffRhs(std::string field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map, FVal exact_value, FVal exact_dot, FVal exact_lap)
boost::function< double(const double, const double, const double)> FVal
std::map< int, BlockData > setOfBlock
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpComputeSlowValue(std::string mass_field, boost::shared_ptr< PreviousData > &data1, boost::shared_ptr< PreviousData > &data2, boost::shared_ptr< PreviousData > &data3, std::map< int, BlockData > &block_map)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData2
boost::shared_ptr< PreviousData > commonData3
boost::shared_ptr< PreviousData > commonData1
boost::shared_ptr< PreviousData > prevData
std::map< int, BlockData > setOfBlock
OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad, boost::shared_ptr< PreviousData > &prev_data, std::map< int, BlockData > &block_map, double &err)
boost::function< double(const double, const double, const double)> FVal
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
VectorDouble dot_values
Rate of values of field "u" at integration points.
MatrixDouble invJac
Inverse of element jacobian.
VectorDouble values
Values of field "u" at integration points.
MatrixDouble grads
Gradients of field "u" at integration points.