23 #ifndef _MIX_TRANPORT_ELEMENT_HPP_
24 #define _MIX_TRANPORT_ELEMENT_HPP_
104 std::vector<int> ids;
108 ids.empty() ? PETSC_NULL : &ids[0],
109 PETSC_COPY_VALUES, &is_local);
110 CHKERR ISAllGather(is_local, is);
111 CHKERR ISDestroy(&is_local);
125 const double y,
const double z,
142 const double y,
const double z,
146 for (
int dd = 0;
dd < 3;
dd++) {
162 const double x,
const double y,
163 const double z,
double &value) {
179 const double y,
const double z,
194 std::map<int, BlockData>
224 const std::string &fluxes_name,
const std::string &values_name) {
259 Mat_Thermal temp_data;
260 CHKERR it->getAttributeDataStructure(temp_data);
262 temp_data.data.Conductivity;
263 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
265 it->meshset, MBTET,
setOfBlocks[it->getMeshsetId()].tEts,
true);
267 setOfBlocks[it->getMeshsetId()].tEts, MBTET,
"MIX");
271 setOfBlocks[it->getMeshsetId()].tEts, 2,
false, skeleton,
272 moab::Interface::UNION);
322 ref_tets = subtract(ref_tets, done_tets);
336 ref_faces = subtract(ref_faces, done_faces);
355 ProblemsManager *prb_mng_ptr;
357 CHKERR prb_mng_ptr->buildProblem(
"MIX",
true);
360 CHKERR prb_mng_ptr->partitionProblem(
"MIX");
361 CHKERR prb_mng_ptr->partitionFiniteElements(
"MIX");
363 CHKERR prb_mng_ptr->partitionGhostDofs(
"MIX");
372 std::vector<EntityHandle> &map_gauss_pts)
385 double *error_flux_ptr;
387 th_error_flux, &fe_ent, 1, (
const void **)&error_flux_ptr);
392 double *error_div_ptr;
394 th_error_div, &fe_ent, 1, (
const void **)&error_div_ptr);
399 double *error_jump_ptr;
401 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
407 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
408 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
409 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
417 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
418 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
419 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
427 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
428 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
429 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
465 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"MIX", &
Aij);
480 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
481 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
483 CHKERR VecGhostUpdateBegin(
D0, INSERT_VALUES, SCATTER_FORWARD);
484 CHKERR VecGhostUpdateEnd(
D0, INSERT_VALUES, SCATTER_FORWARD);
486 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
487 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
490 "MIX",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
498 feTri.getOpPtrVector().clear();
499 feTri.getOpPtrVector().push_back(
500 new OpHOSetContravariantPiolaTransformOnFace3D(
HDIV));
504 CHKERR VecGhostUpdateBegin(
D0, INSERT_VALUES, SCATTER_REVERSE);
505 CHKERR VecGhostUpdateEnd(
D0, INSERT_VALUES, SCATTER_REVERSE);
524 feTri.getOpPtrVector().clear();
525 feTri.getOpPtrVector().push_back(
526 new OpHOSetContravariantPiolaTransformOnFace3D(
HDIV));
531 CHKERR MatAssemblyBegin(
Aij, MAT_FINAL_ASSEMBLY);
532 CHKERR MatAssemblyEnd(
Aij, MAT_FINAL_ASSEMBLY);
533 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
534 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
540 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
541 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_F = %6.4e\n", nrm2_F);
552 CHKERR MatZeroRowsColumnsIS(
Aij, essential_bc_ids, 1,
D0,
F);
553 CHKERR ISDestroy(&essential_bc_ids);
563 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
564 PetscPrintf(PETSC_COMM_WORLD,
"With BC nrm2_F = %6.4e\n", nrm2_F);
577 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
579 CHKERR KSPSetFromOptions(solver);
582 CHKERR KSPDestroy(&solver);
586 CHKERR VecNorm(
D, NORM_2, &nrm2_D);
588 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_D = %6.4e\n", nrm2_D);
590 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
591 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
595 "MIX",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
604 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
605 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
621 feTri.getOpPtrVector().clear();
624 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
625 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
631 std::vector<int> ids;
633 std::vector<double> vals(ids.size(), 0);
641 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
642 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_F = %6.4e\n", nrm2_F);
643 const double eps = 1e-8;
664 feTri.getOpPtrVector().clear();
676 const Problem *problem_ptr;
679 "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
680 "jump^2 = %6.4e error tot^2 = %6.4e\n",
714 const std::string flux_name, Mat aij,
Vec f)
716 flux_name, flux_name,
741 if (
Aij == PETSC_NULL)
743 if (row_data.getIndices().size() == 0)
745 if (col_data.getIndices().size() == 0)
748 int nb_row = row_data.getIndices().size();
749 int nb_col = col_data.getIndices().size();
750 NN.resize(nb_row, nb_col,
false);
754 invK.resize(3, 3,
false);
757 &
invK(0, 0), &
invK(0, 1), &
invK(0, 2), &
invK(1, 0), &
invK(1, 1),
760 auto t_n_hdiv_row = row_data.getFTensor1N<3>();
762 int nb_gauss_pts = row_data.getN().size1();
763 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
771 for (
int kk = 0; kk != nb_row; kk++) {
773 &col_data.getVectorN<3>(gg)(0,
HVEC0),
774 &col_data.getVectorN<3>(gg)(0,
HVEC1),
775 &col_data.getVectorN<3>(gg)(0,
HVEC2), 3);
776 t_row(
j) =
w * t_n_hdiv_row(
i) * t_inv_k(
i,
j);
777 for (
int ll = 0; ll != nb_col; ll++) {
778 NN(kk, ll) += t_row(
j) * t_n_hdiv_col(
j);
786 &col_data.getIndices()[0], &
NN(0, 0), ADD_VALUES);
788 if (row_side != col_side || row_type != col_type) {
789 transNN.resize(nb_col, nb_row);
792 &row_data.getIndices()[0], &
transNN(0, 0),
811 int nb_row = data.getIndices().size();
820 invK.resize(3, 3,
false);
825 &
invK(0, 0), &
invK(0, 1), &
invK(0, 2), &
invK(1, 0), &
invK(1, 1),
828 auto t_n_hdiv = data.getFTensor1N<3>();
830 int nb_gauss_pts = data.getN().size1();
831 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
837 for (
int ll = 0; ll != nb_row; ll++) {
838 Nf[ll] +=
w * t_n_hdiv(
i) * t_inv_k(
i,
j) * t_flux(
j);
859 string val_name_col,
Vec f)
874 if (data.getFieldData().size() == 0)
876 int nb_row = data.getIndices().size();
879 divVec.resize(data.getN().size2() / 3, 0);
880 if (
divVec.size() != data.getIndices().size()) {
882 "data inconsistency");
884 int nb_gauss_pts = data.getN().size1();
887 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
890 for (; gg < nb_gauss_pts; gg++) {
893 v = t_base_diff_hdiv(
i,
i);
919 string flux_name_col, Mat aij,
Vec f)
921 val_name_row, flux_name_col,
951 if (
Aij == PETSC_NULL)
953 if (row_data.getFieldData().size() == 0)
955 if (col_data.getFieldData().size() == 0)
957 int nb_row = row_data.getFieldData().size();
958 int nb_col = col_data.getFieldData().size();
959 NN.resize(nb_row, nb_col);
961 divVec.resize(nb_col,
false);
964 auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
966 int nb_gauss_pts = row_data.getN().size1();
967 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
970 v = t_base_diff_hdiv(
i,
i);
973 noalias(
NN) +=
w * outer_prod(row_data.getN(gg),
divVec);
976 &col_data.getIndices()[0], &
NN(0, 0), ADD_VALUES);
977 transNN.resize(nb_col, nb_row);
980 &row_data.getIndices()[0], &
transNN(0, 0),
988 if (data.getIndices().size() == 0)
990 if (data.getIndices().size() != data.getN().size2()) {
992 "data inconsistency");
994 int nb_row = data.getIndices().size();
997 int nb_gauss_pts = data.getN().size1();
998 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1025 if (data.getFieldData().size() == 0)
1028 int nb_row = data.getFieldData().size();
1029 Nf.resize(nb_row,
false);
1031 int nb_gauss_pts = data.getN().size1();
1032 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1039 noalias(
Nf) +=
w * data.getN(gg) *
flux;
1082 if (data.getFieldData().size() == 0)
1084 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1085 nF.resize(data.getIndices().size());
1087 int nb_gauss_pts = data.getN().size1();
1088 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1089 const double x = getCoordsAtGaussPts()(gg, 0);
1090 const double y = getCoordsAtGaussPts()(gg, 1);
1091 const double z = getCoordsAtGaussPts()(gg, 2);
1096 double w = getGaussPts()(2, gg) * 0.5;
1097 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1099 w * prod(data.getVectorN<3>(gg), getNormalsAtGaussPts(gg)) *
1102 noalias(
nF) +=
w * prod(data.getVectorN<3>(gg), getNormal()) * value;
1105 Vec f = (
F != PETSC_NULL) ?
F : getFEMethod()->ts_F;
1107 &
nF[0], ADD_VALUES);
1141 if (data.getFieldData().size() == 0)
1143 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1144 int nb_dofs = data.getFieldData().size();
1145 int nb_gauss_pts = data.getN().size1();
1146 if (3 * nb_dofs !=
static_cast<int>(data.getN().size2())) {
1148 "wrong number of dofs");
1150 NN.resize(nb_dofs, nb_dofs);
1159 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1161 normal_ptr = &getNormalsAtGaussPts(0)[0];
1164 normal_ptr = &getNormal()[0];
1170 auto t_n_hdiv_row = data.getFTensor1N<3>();
1175 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1178 const double x = getCoordsAtGaussPts()(gg, 0);
1179 const double y = getCoordsAtGaussPts()(gg, 1);
1180 const double z = getCoordsAtGaussPts()(gg, 2);
1187 double w = getGaussPts()(2, gg);
1189 nrm2 = sqrt(t_normal(
i) * t_normal(
i));
1194 for (
int ll = 0; ll != nb_dofs; ll++) {
1197 &data.getVectorN<3>(gg)(0,
HVEC0),
1198 &data.getVectorN<3>(gg)(0,
HVEC1),
1199 &data.getVectorN<3>(gg)(0,
HVEC2), 3);
1200 for (
int kk = 0; kk <= ll; kk++) {
1201 NN(ll, kk) +=
w * t_n_hdiv_row(
i) * t_n_hdiv_col(
i);
1205 Nf[ll] +=
w * t_n_hdiv_row(
i) * t_normal(
i) *
flux / nrm2;
1210 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1212 nrm2 = sqrt(t_normal(
i) * t_normal(
i));
1216 cTx.
bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
1223 CHKERR VecSetOption(
cTx.
D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1225 &*
Nf.begin(), INSERT_VALUES);
1226 for (
int dd = 0;
dd != nb_dofs; ++
dd)
1227 data.getFieldDofs()[
dd]->getFieldData() =
Nf[
dd];
1242 const std::string val_name =
"VALUES")
1252 if (data.getFieldData().size() == 0)
1255 int nb_gauss_pts = data.getN().size1();
1257 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1259 inner_prod(trans(data.getN(gg)), data.getFieldData());
1275 const std::string val_name =
"VALUES")
1284 if (data.getFieldData().size() == 0)
1286 int nb_gauss_pts = data.getDiffN().size1();
1288 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1289 ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1291 noalias(values_grad_at_gauss_pts) =
1292 prod(trans(data.getDiffN(gg)), data.getFieldData());
1307 const std::string field_name)
1317 if (data.getFieldData().size() == 0)
1319 int nb_gauss_pts = data.getDiffN().size1();
1320 int nb_dofs = data.getFieldData().size();
1323 if (
type == MBTRI && side == 0) {
1330 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
1332 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1334 v = t_base_diff_hdiv(
i,
i);
1338 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1341 prod(trans(data.getVectorN<3>(gg)), data.getFieldData());
1373 invK.resize(3, 3,
false);
1374 int nb_gauss_pts = data.getN().size1();
1377 Tag th_error_flux, th_error_div;
1379 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
1380 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1381 double *error_flux_ptr;
1383 th_error_flux, &fe_ent, 1, (
const void **)&error_flux_ptr);
1386 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
1387 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1388 double *error_div_ptr;
1390 th_error_div, &fe_ent, 1, (
const void **)&error_div_ptr);
1394 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1395 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1396 double *error_jump_ptr;
1398 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
1399 *error_jump_ptr = 0;
1402 const double h = pow(
getVolume() * 12 / std::sqrt(2), (
double)1 / 3);
1404 for (
int ff = 0; ff != 4; ff++) {
1407 double *error_face_jump_ptr;
1409 th_error_jump, &face, 1, (
const void **)&error_face_jump_ptr);
1410 *error_face_jump_ptr = (1 / sqrt(
h)) * sqrt(*error_face_jump_ptr);
1411 *error_face_jump_ptr = pow(*error_face_jump_ptr, 2);
1412 *error_jump_ptr += *error_face_jump_ptr;
1415 *error_flux_ptr = 0;
1418 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1429 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1431 ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1434 prod(
invK, flux_at_gauss_pt) + values_grad_at_gauss_pts;
1437 *error_div_ptr =
h * sqrt(*error_div_ptr);
1438 *error_div_ptr = pow(*error_div_ptr, 2);
1439 cTx.
errorMap[sqrt(*error_flux_ptr + *error_div_ptr + *error_jump_ptr)] =
1480 if (data.getFieldData().size() == 0)
1482 int nb_gauss_pts = data.getN().size1();
1483 valMap[getFaceSense()].resize(nb_gauss_pts);
1484 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1485 valMap[getFaceSense()][gg] =
1486 inner_prod(trans(data.getN(gg)), data.getFieldData());
1505 if (
type == MBTRI) {
1506 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1511 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1512 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1513 double *error_jump_ptr;
1515 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
1516 *error_jump_ptr = 0;
1532 int nb_gauss_pts = data.getN().size1();
1535 if (
valMap.size() == 1) {
1536 if (
static_cast<int>(
valMap.begin()->second.size()) != nb_gauss_pts) {
1538 "wrong number of integration points");
1540 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
1541 const double x = getCoordsAtGaussPts()(gg, 0);
1542 const double y = getCoordsAtGaussPts()(gg, 1);
1543 const double z = getCoordsAtGaussPts()(gg, 2);
1546 double w = getGaussPts()(2, gg);
1547 if (
static_cast<int>(getNormalsAtGaussPts().size1()) ==
1549 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1553 *error_jump_ptr +=
w * pow(value -
valMap.begin()->second[gg], 2);
1555 }
else if (
valMap.size() == 2) {
1556 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
1557 double w = getGaussPts()(2, gg);
1558 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1559 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1564 *error_jump_ptr +=
w * pow(
delta, 2);
1568 "data inconsistency, wrong number of neighbors "
1569 "valMap.size() = %d",
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
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.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define BITREFLEVEL_SIZE
max number of refinements
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MAT_THERMALSET
block name is "MAT_THERMAL"
@ 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 ...
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual EntityHandle get_finite_element_meshset(const std::string &name) const =0
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreOnSideSwitch< 0 > VolumeElementForcesAndSourcesCoreOnSide
Volume element used to integrate on skeleton.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
MatrixBoundedArray< double, 9 > MatrixDouble3by3
implementation of Data Operators for Forces and Sources
DeprecatedCoreInterface Interface
constexpr auto VecSetValues
constexpr auto MatSetValues
double flux
impulse magnitude
double h
convective heat coefficient
static constexpr double delta
data for evaluation of het conductivity and heat capacity elements
Range tEts
constatins elements in block set
definition of surface element
MyTriFE(MoFEM::Interface &m_field)
definition of volume element
MyVolumeFE(MoFEM::Interface &m_field)
MixTransportElement & cTx
OpDivTauU_HdivL2(MixTransportElement &ctx, const std::string flux_name_row, string val_name_col, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate error evaluator
OpError(MixTransportElement &ctx, const std::string field_name)
MixTransportElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Evaluate boundary conditions on fluxes.
FTensor::Index< 'i', 3 > i
MixTransportElement & cTx
OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate flux at integration point
OpFluxDivergenceAtGaussPts(MixTransportElement &ctx, const std::string field_name)
MixTransportElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate source therms, i.e. .
OpL2Source(MixTransportElement &ctx, const std::string val_name, Vec f)
MixTransportElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpPostProc(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
MixTransportElement & cTx
OpRhsBcOnValues(MixTransportElement &ctx, const std::string fluxes_name, Vec f)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
calculate values on adjacent tetrahedra to face
map< int, VectorDouble > & valMap
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpVolSide(map< int, VectorDouble > &val_map)
calculate jump on entities
OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
map< int, VectorDouble > valMap
VolumeElementForcesAndSourcesCoreOnSide volSideFe
volume element to get values from adjacent tets to face
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MixTransportElement & cTx
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
MixTransportElement & cTx
OpTauDotSigma_HdivHdiv(MixTransportElement &ctx, const std::string flux_name, Mat aij, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Assemble matrix.
MixTransportElement & cTx
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
OpVDivSigma_L2Hdiv(MixTransportElement &ctx, const std::string val_name_row, string flux_name_col, Mat aij, Vec f)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate values at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpValuesAtGaussPts(MixTransportElement &ctx, const std::string val_name="VALUES")
MixTransportElement & cTx
virtual ~OpValuesAtGaussPts()
Calculate gradients of values at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpValuesGradientAtGaussPts(MixTransportElement &ctx, const std::string val_name="VALUES")
MixTransportElement & cTx
virtual ~OpValuesGradientAtGaussPts()
MoFEMErrorCode postProc(const string out_file)
Post process results.
MoFEMErrorCode destroyMatrices()
destroy matrices
MoFEMErrorCode evaluateError()
Calculate error on elements.
MoFEMErrorCode solveLinearProblem()
solve problem
MoFEMErrorCode getDirichletBCIndices(IS *is)
get dof indices where essential boundary conditions are applied
MoFEMErrorCode createMatrices()
create matrices
MoFEM::Interface & mField
MatrixDouble valuesGradientAtGaussPts
gradients at integration points on element
virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode calculateResidual()
calculate residual
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
VectorDouble valuesAtGaussPts
values at integration points on element
map< double, EntityHandle > errorMap
MyTriFE feTri
Instance of surface element.
virtual MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
evaluate natural (Dirichlet) boundary conditions
MyVolumeFE feVol
Instance of volume element.
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
Add fields to database.
virtual MoFEMErrorCode getResistivity(const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
natural (Dirichlet) boundary conditions (set values)
MixTransportElement(MoFEM::Interface &m_field)
construction of this data structure
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
VectorDouble divergenceAtGaussPts
divergence at integration points on element
virtual ~MixTransportElement()
destructor
virtual MoFEMErrorCode getSource(const EntityHandle ent, const double x, const double y, const double z, double &flux)
set source term
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
MoFEMErrorCode buildProblem(BitRefLevel &ref_level)
Build problem.
virtual int get_comm_size() const =0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Mat & ts_B
Preconditioner for ts_A.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
double getVolume() const
element volume (linear geometry)
VolumeElementForcesAndSourcesCoreBase * getVolumeFE() const
return pointer to Generic Volume Finite Element object
Volume finite element with switches.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.