1#ifndef __RDOPERATORS_HPP__
2#define __RDOPERATORS_HPP__
9using FaceEle = MoFEM::FaceElementForcesAndSourcesCoreSwitch<0>;
11using BoundaryEle = MoFEM::EdgeElementForcesAndSourcesCoreSwitch<0>;
16using EntData = DataForcesAndSourcesCore::EntData;
44 jac.resize(2, 2,
false);
82 boost::shared_ptr<PreviousData> &data1,
83 boost::shared_ptr<PreviousData> &data2,
84 boost::shared_ptr<PreviousData> &data3,
85 std::map<int, BlockData> &block_map)
91 boost::shared_ptr<VectorDouble> slow_value_ptr1(
93 boost::shared_ptr<VectorDouble> slow_value_ptr2(
95 boost::shared_ptr<VectorDouble> slow_value_ptr3(
101 const int nb_integration_pts =
getGaussPts().size2();
102 if (type == MBVERTEX) {
103 vec1.resize(nb_integration_pts,
false);
104 vec2.resize(nb_integration_pts,
false);
105 vec3.resize(nb_integration_pts,
false);
110 const int nb_dofs = data.getIndices().size();
113 auto find_block_data = [&]() {
117 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
118 block_raw_ptr = &
m.second;
122 return block_raw_ptr;
125 auto block_data_ptr = find_block_data();
129 auto &block_data = *block_data_ptr;
131 const int nb_integration_pts =
getGaussPts().size2();
141 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
142 t_slow_values1 = block_data.r1 * t_mass_values1 *
143 (1.0 - block_data.a11 * t_mass_values1 -
144 block_data.a12 * t_mass_values2 -
145 block_data.a13 * t_mass_values3);
146 t_slow_values2 = block_data.r2 * t_mass_values2 *
147 (1.0 - block_data.a21 * t_mass_values1 -
148 block_data.a22 * t_mass_values2 -
149 block_data.a23 * t_mass_values3);
151 t_slow_values3 = block_data.r3 * t_mass_values3 *
152 (1.0 - block_data.a31 * t_mass_values1 -
153 block_data.a32 * t_mass_values2 -
154 block_data.a33 * t_mass_values3);
182 int nb_dofs = data.getIndices().size();
189 int size2 = data.getN().size2();
190 if (3 * nb_dofs !=
static_cast<int>(data.getN().size2()))
192 "wrong number of dofs");
193 nN.resize(nb_dofs, nb_dofs,
false);
194 nF.resize(nb_dofs,
false);
198 auto t_row_tau = data.getFTensor1N<3>();
200 auto dir = getDirection();
201 double len = sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
209 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
210 const double a = t_w * vol;
211 for (
int rr = 0; rr != nb_dofs; rr++) {
212 auto t_col_tau = data.getFTensor1N<3>(gg, 0);
214 for (
int cc = 0; cc != nb_dofs; cc++) {
215 nN(rr, cc) +=
a * (t_row_tau(
i) * t_normal(
i)) *
216 (t_col_tau(
i) * t_normal(
i));
227 for (
auto &dof : data.getFieldDofs()) {
228 dof->getFieldData() =
nF[dof->getEntDofIdx()];
247 Range &internal_edge_ents)
257 int nb_dofs = data.getIndices().size();
264 int size2 = data.getN().size2();
265 if (3 * nb_dofs !=
static_cast<int>(data.getN().size2()))
267 "wrong number of dofs");
268 vecF.resize(nb_dofs,
false);
272 auto t_row_v_base = data.getFTensor0N();
281 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
282 const double a = t_w * vol;
285 for (
int rr = 0; rr != nb_dofs; ++rr) {
287 vecF[rr] +=
a * trace_val * t_row_v_base;
327 int nb_dofs = data.getFieldData().size();
333 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
335 "wrong number of dofs");
336 nN.resize(nb_dofs, nb_dofs,
false);
337 nF.resize(nb_dofs,
false);
341 auto t_row_mass = data.getFTensor0N();
345 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
346 const double a = t_w * vol;
348 for (
int rr = 0; rr != nb_dofs; rr++) {
349 auto t_col_mass = data.getFTensor0N(gg, 0);
351 for (
int cc = 0; cc != nb_dofs; cc++) {
352 nN(rr, cc) +=
a * t_row_mass * t_col_mass;
363 for (
auto &dof : data.getFieldDofs()) {
364 dof->getFieldData() =
nF[dof->getEntDofIdx()];
381 typedef boost::function<
double(
const double,
const double,
const double)>
383 typedef boost::function<FTensor::Tensor1<double, 3>(
384 const double,
const double,
const double)>
387 boost::shared_ptr<PreviousData> &common_data,
402 const int nb_dofs = data.getIndices().size();
405 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
407 "wrong number of dofs");
408 vecF.resize(nb_dofs,
false);
409 mat.resize(nb_dofs, nb_dofs,
false);
412 const int nb_integration_pts =
getGaussPts().size2();
415 auto t_row_v_base = data.getFTensor0N();
421 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
422 const double a = vol * t_w;
424 double u_dot =
exactDot(t_coords(
NX), t_coords(
NY), ct);
425 double u_lap =
exactLap(t_coords(
NX), t_coords(
NY), ct);
429 for (
int rr = 0; rr != nb_dofs; ++rr) {
430 auto t_col_v_base = data.getFTensor0N(gg, 0);
431 vecF[rr] +=
a * t_slow_value * t_row_v_base;
433 for (
int cc = 0; cc != nb_dofs; ++cc) {
434 mat(rr, cc) +=
a * t_row_v_base * t_col_v_base;
477 const int nb_dofs = data.getIndices().size();
480 EntityHandle row_side_ent = data.getFieldDofs()[0]->getEnt();
486 vecF.resize(nb_dofs,
false);
488 const int nb_integration_pts =
getGaussPts().size2();
489 auto t_tau_base = data.getFTensor1N<3>();
491 auto dir = getDirection();
496 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
497 const double a = t_w;
498 for (
int rr = 0; rr != nb_dofs; ++rr) {
525 boost::shared_ptr<PreviousData> &data,
526 std::map<int, BlockData> &block_map)
533 const int nb_dofs = data.getIndices().size();
535 auto find_block_data = [&]() {
539 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
540 block_raw_ptr = &
m.second;
544 return block_raw_ptr;
547 auto block_data_ptr = find_block_data();
550 auto &block_data = *block_data_ptr;
552 vecF.resize(nb_dofs,
false);
555 const int nb_integration_pts =
getGaussPts().size2();
556 auto t_flux_value = getFTensor1FromMat<3>(
commonData->flux_values);
558 auto t_tau_base = data.getFTensor1N<3>();
560 auto t_tau_grad = data.getFTensor2DiffN<3, 2>();
565 for (
int gg = 0; gg < nb_integration_pts; ++gg) {
567 const double K =
B_epsilon + (block_data.B0 +
B * t_mass_value);
568 const double K_inv = 1. /
K;
569 const double a = vol * t_w;
570 for (
int rr = 0; rr < nb_dofs; ++rr) {
571 double div_base = t_tau_grad(0, 0) + t_tau_grad(1, 1);
572 vecF[rr] += (K_inv * t_tau_base(
i) * t_flux_value(
i) -
573 div_base * t_mass_value) *
599 typedef boost::function<
double(
const double,
const double,
const double)>
602 boost::shared_ptr<PreviousData> &data,
603 std::map<int, BlockData> &block_map,
617 const int nb_dofs = data.getIndices().size();
620 auto find_block_data = [&]() {
624 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
625 block_raw_ptr = &
m.second;
629 return block_raw_ptr;
632 auto block_data_ptr = find_block_data();
635 auto &block_data = *block_data_ptr;
638 vecF.resize(nb_dofs,
false);
640 const int nb_integration_pts =
getGaussPts().size2();
643 auto t_row_v_base = data.getFTensor0N();
649 for (
int gg = 0; gg < nb_integration_pts; ++gg) {
650 const double a = vol * t_w;
651 double u_dot =
exactDot(t_coords(
NX), t_coords(
NY), ct);
652 double u_lap =
exactLap(t_coords(
NX), t_coords(
NY), ct);
654 double f = u_dot - block_data.B0 * u_lap;
655 for (
int rr = 0; rr < nb_dofs; ++rr) {
656 vecF[rr] += (t_row_v_base * (t_mass_dot + t_flux_div)) *
a;
693 std::map<int, BlockData> &block_map)
704 const int nb_row_dofs = row_data.getIndices().size();
705 const int nb_col_dofs = col_data.getIndices().size();
707 if (nb_row_dofs && nb_col_dofs) {
708 auto find_block_data = [&]() {
712 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
713 block_raw_ptr = &
m.second;
717 return block_raw_ptr;
720 auto block_data_ptr = find_block_data();
723 auto &block_data = *block_data_ptr;
725 mat.resize(nb_row_dofs, nb_col_dofs,
false);
727 const int nb_integration_pts =
getGaussPts().size2();
730 auto t_row_tau_base = row_data.getFTensor1N<3>();
735 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
736 const double a = vol * t_w;
737 const double K =
B_epsilon + (block_data.B0 +
B * t_mass_value);
738 const double K_inv = 1. /
K;
739 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
740 auto t_col_tau_base = col_data.getFTensor1N<3>(gg, 0);
741 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
742 mat(rr, cc) += (K_inv * t_row_tau_base(
i) * t_col_tau_base(
i)) *
a;
752 if (row_side != col_side || row_type != col_type) {
753 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
774 boost::shared_ptr<PreviousData> &data,
775 std::map<int, BlockData> &block_map)
786 const int nb_row_dofs = row_data.getIndices().size();
787 const int nb_col_dofs = col_data.getIndices().size();
789 if (nb_row_dofs && nb_col_dofs) {
790 auto find_block_data = [&]() {
794 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
795 block_raw_ptr = &
m.second;
799 return block_raw_ptr;
802 auto block_data_ptr = find_block_data();
805 auto &block_data = *block_data_ptr;
806 mat.resize(nb_row_dofs, nb_col_dofs,
false);
808 const int nb_integration_pts =
getGaussPts().size2();
810 auto t_row_tau_base = row_data.getFTensor1N<3>();
812 auto t_row_tau_grad = row_data.getFTensor2DiffN<3, 2>();
814 auto t_flux_value = getFTensor1FromMat<3>(
commonData->flux_values);
817 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
818 const double a = vol * t_w;
819 const double K =
B_epsilon + (block_data.B0 +
B * t_mass_value);
820 const double K_inv = 1. /
K;
821 const double K_diff =
B;
823 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
824 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
825 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
826 double div_row_base = t_row_tau_grad(0, 0) + t_row_tau_grad(1, 1);
827 mat(rr, cc) += (-(t_row_tau_base(
i) * t_flux_value(
i) * K_inv *
828 K_inv * K_diff * t_col_v_base) -
829 (div_row_base * t_col_v_base)) *
865 const int nb_row_dofs = row_data.getIndices().size();
866 const int nb_col_dofs = col_data.getIndices().size();
868 if (nb_row_dofs && nb_col_dofs) {
869 mat.resize(nb_row_dofs, nb_col_dofs,
false);
871 const int nb_integration_pts =
getGaussPts().size2();
873 auto t_row_v_base = row_data.getFTensor0N();
876 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
877 const double a = vol * t_w;
878 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
879 auto t_col_tau_grad = col_data.getFTensor2DiffN<3, 2>(gg, 0);
880 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
881 double div_col_base = t_col_tau_grad(0, 0) + t_col_tau_grad(1, 1);
882 mat(rr, cc) += (t_row_v_base * div_col_base) *
a;
913 const int nb_row_dofs = row_data.getIndices().size();
914 const int nb_col_dofs = col_data.getIndices().size();
915 if (nb_row_dofs && nb_col_dofs) {
917 mat.resize(nb_row_dofs, nb_col_dofs,
false);
919 const int nb_integration_pts =
getGaussPts().size2();
921 auto t_row_v_base = row_data.getFTensor0N();
927 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
928 const double a = vol * t_w;
930 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
931 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
932 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
933 mat(rr, cc) += (ts_a * t_row_v_base * t_col_v_base) *
a;
943 if (row_side != col_side || row_type != col_type) {
944 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
958 typedef boost::function<
double(
const double,
const double,
const double)>
960 typedef boost::function<FTensor::Tensor1<double, 3>(
961 const double,
const double,
const double)>
966 boost::shared_ptr<PreviousData> &prev_data,
967 std::map<int, BlockData> &block_map,
979 const int nb_dofs = data.getFieldData().size();
982 auto find_block_data = [&]() {
986 if (
m.second.block_ents.find(fe_ent) !=
m.second.block_ents.end()) {
987 block_raw_ptr = &
m.second;
991 return block_raw_ptr;
994 auto block_data_ptr = find_block_data();
998 auto &block_data = *block_data_ptr;
1002 auto t_flux_value = getFTensor1FromMat<3>(
prevData->flux_values);
1006 data.getFieldData().clear();
1008 const int nb_integration_pts =
getGaussPts().size2();
1017 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1018 const double a = vol * t_w;
1019 double mass_exact =
exactVal(t_coords(
NX), t_coords(
NY), ct);
1021 t_exact_flux(
i) = - block_data.B0 *
exactGrad(t_coords(
NX), t_coords(
NY), ct)(
i);
1022 t_flux_error(0) = t_flux_value(0) - t_exact_flux(0);
1023 t_flux_error(1) = t_flux_value(1) - t_exact_flux(1);
1024 t_flux_error(2) = t_flux_value(2) - t_exact_flux(2);
1025 double local_error = pow(mass_exact - t_mass_value, 2) + t_flux_error(
i) * t_flux_error(
i);
1027 data.getFieldData()[0] +=
a * local_error;
1028 eRror +=
a * local_error;
1038 data.getFieldDofs()[0]->getFieldData() = data.getFieldData()[0];
1068 boost::shared_ptr<PostProcFaceOnRefinedMesh> &post_proc,
double &err)
1084 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
1087 CHKERR VecCreateMPI(
cOmm, 1, PETSC_DECIDE, &error_per_proc);
1088 auto get_global_error = [&]() {
1093 CHKERR get_global_error();
1094 CHKERR VecAssemblyBegin(error_per_proc);
1095 CHKERR VecAssemblyEnd(error_per_proc);
1097 CHKERR VecSum(error_per_proc, &error_sum);
1098 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Error : %3.4e \n",
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
#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.
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
double compute_natu_bc(const double x, const double y, const double z)
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< 0 > BoundaryEle
double compute_init_val(const double x, const double y, const double z)
MoFEM::FaceElementForcesAndSourcesCoreSwitch< 0 > FaceEle
DataForcesAndSourcesCore::EntData EntData
double compute_essen_bc(const double x, const double y, const double z)
bool sYmm
If true assume that matrix is symmetric structure.
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.
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 postProcess()
function is run at the end of loop
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)
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProc
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsTauTau(std::string flux_field, boost::shared_ptr< PreviousData > &commonData, std::map< int, BlockData > &block_map)
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)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
OpAssembleLhsTauV(std::string flux_field, std::string mass_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsVV(std::string mass_field)
OpAssembleNaturalBCRhsTau(std::string flux_field, Range &natural_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonData
boost::function< double(const double, const double, const double)> FVal
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &common_data, FVal exact_value, FVal exact_dot, FVal exact_lap)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
boost::shared_ptr< PreviousData > commonData
OpAssembleStiffRhsV(std::string flux_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map, FVal exact_value, FVal exact_dot, FVal exact_lap)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
std::map< int, BlockData > setOfBlock
boost::function< double(const double, const double, const double)> FVal
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonData2
boost::shared_ptr< PreviousData > commonData1
std::map< int, BlockData > setOfBlock
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)
boost::shared_ptr< PreviousData > commonData3
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > prevData
std::map< int, BlockData > setOfBlock
boost::function< double(const double, const double, const double)> FVal
OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad, boost::shared_ptr< PreviousData > &prev_data, std::map< int, BlockData > &block_map, double &err)
Range & essential_bd_ents
OpEssentialBC(const std::string &flux_field, Range &essential_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpInitialMass(const std::string &mass_field, Range &inner_surface)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< double(const double, const double, const double)> ExactFunVal
boost::function< double(const double)> FVal
OpSkeletonSource(const std::string &mass_field, FVal skeleton_fun, ExactFunVal smooth_fun, Range &internal_edge_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)