9using namespace boost::numeric;
14 MoFEMErrorCode(
const boost::shared_ptr<MoFEM::NumeredDofEntity> &dof)>
20 for (
auto eit = ents.pair_begin(); eit != ents.pair_end(); ++eit) {
22 auto lo_dit = dofs_by_uid.lower_bound(
24 auto hi_dit = dofs_by_uid.upper_bound(
27 for (; lo_dit != hi_dit; ++lo_dit) {
29 if (dof->getHasLocalIndex())
39 Mat Aij, Vec X, Vec
F,
41 : mField(m_field), fieldName(
field_name), blocksetName(blockset_name),
54 : mField(m_field), fieldName(
field_name), blocksetName(blockset_name),
65 std::vector<DataFromBc> &bc_data) {
74 CHKERR bc_data.back().getBcData(mydata, &(*it));
81 std::vector<double> mydata;
82 CHKERR bc_data.back().getBcData(mydata, &(*it));
100 const double angle = std::sqrt(t_omega(
i) * t_omega(
i));
101 if (std::abs(angle) < 1e-18)
105 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
106 const double a = sin(angle) / angle;
107 const double ss_2 = sin(angle / 2.);
108 const double b = 2. * ss_2 * ss_2 / (angle * angle);
109 t_R(
i,
j) +=
a * t_Omega(
i,
j);
110 t_R(
i,
j) += b * t_Omega(
i,
k) * t_Omega(
k,
j);
122 FTensor1 t_coords(coords[0], coords[1], coords[2]);
123 const double a = std::sqrt(t_normal(
i) * t_normal(
i));
124 t_omega(
i) = t_normal(
i) * (theta /
a);
127 t_delta(
i) = t_centr(
i) - t_coords(
i);
128 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
134 std::vector<DataFromBc> &bc_data) {
139 if (it->getName().compare(0, 8,
"ROTATION") == 0) {
141 std::vector<double> mydata;
142 CHKERR bc_data.back().getEntitiesFromBc(
mField, &(*it));
143 CHKERR it->getAttributes(mydata);
144 if (mydata.size() < 6) {
146 "6 attributes are required for Rotation (3 center coords + 3 "
147 "angles, (+ 3 optional) flags for xyz)");
149 for (
int ii : {0, 1, 2}) {
150 bc_data.back().bc_flags[ii] = 1;
151 bc_data.back().t_centr(ii) = mydata[ii + 0];
152 bc_data.back().t_normal(ii) = mydata[ii + 3];
154 if (mydata.size() > 8)
155 for (
int ii : {0, 1, 2})
156 bc_data.back().bc_flags[ii] = mydata[6 + ii];
158 bc_data.back().scaled_values[0] = bc_data.back().initial_values[0] =
159 bc_data.back().t_normal.l2();
160 bc_data.back().is_rotation =
true;
202 std::vector<DataFromBc> bcData;
206 for (
auto &bc_it : bcData) {
212 auto for_each_dof = [&](
auto &dof) {
214 if (dof->getEntType() == MBVERTEX) {
216 if (dof->getDofCoeffIdx() == 0 && bc_it.bc_flags[0]) {
218 mapZeroRows[dof->getPetscGlobalDofIdx()] = bc_it.scaled_values[0];
219 dof->getFieldData() =
mapZeroRows[dof->getPetscGlobalDofIdx()];
221 if (dof->getDofCoeffIdx() == 1 && bc_it.bc_flags[1]) {
222 mapZeroRows[dof->getPetscGlobalDofIdx()] = bc_it.scaled_values[1];
223 dof->getFieldData() =
mapZeroRows[dof->getPetscGlobalDofIdx()];
225 if (dof->getDofCoeffIdx() == 2 && bc_it.bc_flags[2]) {
226 mapZeroRows[dof->getPetscGlobalDofIdx()] = bc_it.scaled_values[2];
227 dof->getFieldData() =
mapZeroRows[dof->getPetscGlobalDofIdx()];
231 if (dof->getDofCoeffIdx() == 0 && bc_it.bc_flags[0]) {
233 dof->getFieldData() = 0;
235 if (dof->getDofCoeffIdx() == 1 && bc_it.bc_flags[1]) {
237 dof->getFieldData() = 0;
239 if (dof->getDofCoeffIdx() == 2 && bc_it.bc_flags[2]) {
241 dof->getFieldData() = 0;
249 bc_it.bc_ents[
dim], for_each_dof);
291 dIag, PETSC_NULL, PETSC_NULL);
296 for (std::vector<int>::iterator vit =
dofsIndices.begin();
311 const double *a_snes_x;
313 auto &dofs_by_glob_idx =
317 auto dof_it = dofs_by_glob_idx.find(git);
318 if (dof_it != dofs_by_glob_idx.end()) {
319 dofsXValues[idx] = a_snes_x[dof_it->get()->getPetscLocalDofIdx()];
323 "Dof with global %d id not found", git);
332 for (std::vector<int>::iterator vit =
dofsIndices.begin();
358 dIag, PETSC_NULL, PETSC_NULL);
373 std::vector<DataFromBc> bcData;
382 for (
auto &bc_it : bcData) {
388 auto for_each_dof = [&](
auto &dof) {
394 if (!dof->getDofCoeffIdx()) {
398 if (eit != field_ents_by_uid.end())
399 noalias(coords) = (*eit)->getEntFieldData();
402 &*coords.data().begin());
405 if (dof->getDofCoeffIdx() == 0 && bc_it.bc_flags[0]) {
407 coords[0] + bc_it.scaled_values[0];
409 dof->getFieldData() =
mapZeroRows[dof->getPetscGlobalDofIdx()];
411 if (dof->getDofCoeffIdx() == 1 && bc_it.bc_flags[1]) {
413 coords[1] + bc_it.scaled_values[1];
414 dof->getFieldData() =
mapZeroRows[dof->getPetscGlobalDofIdx()];
416 if (dof->getDofCoeffIdx() == 2 && bc_it.bc_flags[2]) {
418 coords[2] + bc_it.scaled_values[2];
419 dof->getFieldData() =
mapZeroRows[dof->getPetscGlobalDofIdx()];
422 if (dof->getDofCoeffIdx() == 0 && bc_it.bc_flags[0]) {
424 dof->getFieldData() = 0;
426 if (dof->getDofCoeffIdx() == 1 && bc_it.bc_flags[1]) {
428 dof->getFieldData() = 0;
430 if (dof->getDofCoeffIdx() == 2 && bc_it.bc_flags[2]) {
432 dof->getFieldData() = 0;
441 bc_it.bc_ents[
dim], for_each_dof);
443 auto fix_field_dof = [&](
auto &dof) {
452 bc_it.bc_ents[
dim], fix_field_dof);
461 for (
int ii = 0; mit !=
mapZeroRows.end(); mit++, ii++) {
473 auto insert_temp_bc = [&](
double &
temp,
auto &it) {
485 moab::Interface::UNION);
486 ents.insert(_edges.begin(), _edges.end());
491 ents.insert(_nodes.begin(), _nodes.end());
494 auto for_each_dof = [&](
auto &dof) {
497 if (dof->getEntType() == MBVERTEX) {
498 mapZeroRows[dof->getPetscGlobalDofIdx()] = temp_2[0];
515 bool is_blockset_defined =
false;
519 std::vector<double> mydata;
520 CHKERR it->getAttributes(mydata);
524 "missing temperature attribute");
526 double my_temp = mydata[0];
527 CHKERR insert_temp_bc(my_temp, it);
529 is_blockset_defined =
true;
535 if (!is_blockset_defined) {
539 CHKERR it->getBcDataStructure(mydata);
541 scaled_values[0] = mydata.
data.value1;
542 CHKERR insert_temp_bc(scaled_values[0], it);
548 std::map<DofIdx, FieldData>::iterator mit =
mapZeroRows.begin();
564 auto for_each_dof = [&](
auto &dof) {
603 dIag, PETSC_NULL, PETSC_NULL);
609 for (std::vector<int>::iterator vit =
dofsIndices.begin();
625 for (std::vector<int>::iterator vit =
dofsIndices.begin();
638 dIag, PETSC_NULL, PETSC_NULL);
653 const int nb_coefficients = field_ptr->getNbOfCoeffs();
655 bcDataPtr = boost::make_shared<vector<DataFromBc>>();
659 auto get_dim = [&](
const Range &ents) {
660 for (
auto d : {3, 2, 1})
661 if (ents.num_of_dimension(d))
666 auto get_adj_ents = [&](
const Range &ents) {
669 const auto dim = get_dim(ents);
670 for (
size_t d = 1; d <
dim; ++d)
672 moab::Interface::UNION);
678 auto remove_dofs_from_dirichlet_data = [&]() {
681 for (
auto &ents : bc_it.bc_ents)
682 for (
int i = 0;
i != nb_coefficients;
i++)
683 if (bc_it.bc_flags[
i])
691 auto remove_dofs_from_dirichlet_data_non_distributed = [&]() {
694 for (
auto &ents : bc_it.bc_ents)
695 for (
int i = 0;
i != nb_coefficients;
i++)
696 if (bc_it.bc_flags[
i])
697 CHKERR prb_mng->removeDofsOnEntitiesNotDistributed(
705 CHKERR remove_dofs_from_dirichlet_data();
707 CHKERR remove_dofs_from_dirichlet_data_non_distributed();
721 for (
auto &ents : bc_it.bc_ents)
722 all_bc_ents.merge(ents);
739 this->bc_flags[1] = (
int)mydata.
data.flag2;
748 if (mydata.size() < 6) {
750 "six attributes are required for given BC blockset (3 values + "
753 for (
unsigned int ii = 0; ii < 3; ii++) {
755 this->
bc_flags[ii] = (int)mydata[ii + 3];
770 moab::Interface::UNION);
771 bc_ents[1].insert(edges.begin(), edges.end());
776 bc_ents[0].insert(nodes.begin(), nodes.end());
789 CHKERR VecGetArrayRead(internal, &array);
794 std::vector<int> ghosts(nb_coefficients);
795 for (
int g = 0;
g != nb_coefficients; ++
g)
802 &*ghosts.begin(), &
v);
807 const int id = it->getMeshsetId();
809 reaction_vec.resize(nb_coefficients);
810 reaction_vec.clear();
819 verts.insert(nodes.begin(), nodes.end());
822 auto for_each_dof = [&](
auto &dof) {
824 reaction_vec[dof->getDofCoeffIdx()] += array[dof->getPetscLocalDofIdx()];
830 verts, for_each_dof);
834 CHKERR VecGetArray(
v, &res_array);
835 for (
int dd = 0; dd != reaction_vec.size(); ++dd)
836 res_array[dd] = reaction_vec[dd];
837 CHKERR VecRestoreArray(
v, &res_array);
839 CHKERR VecGetArray(
v, &res_array);
840 for (
int dd = 0; dd != reaction_vec.size(); ++dd)
841 reaction_vec[dd] = res_array[dd];
842 CHKERR VecRestoreArray(
v, &res_array);
845 CHKERR VecGhostUpdateBegin(
v, ADD_VALUES, SCATTER_REVERSE);
846 CHKERR VecGhostUpdateEnd(
v, ADD_VALUES, SCATTER_REVERSE);
847 CHKERR VecGhostUpdateBegin(
v, INSERT_VALUES, SCATTER_FORWARD);
848 CHKERR VecGhostUpdateEnd(
v, INSERT_VALUES, SCATTER_FORWARD);
851 CHKERR VecRestoreArrayRead(internal, &array);
auto get_rotation_from_vector(FTensor1 &t_omega)
static MoFEMErrorCode set_numered_dofs_on_ents(const Problem *problem_ptr, const FieldBitNumber bit_number, Range &ents, boost::function< MoFEMErrorCode(const boost::shared_ptr< MoFEM::NumeredDofEntity > &dof)> for_each_dof)
auto get_displacement(VectorDouble3 &coords, FTensor1 t_centr, FTensor1 t_normal, double theta)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
#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_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 ...
void temp(int x, int y=10)
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, member< FieldEntity, UId, &FieldEntity::localUId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FieldEntity::interface_type_RefEntity, EntityHandle, &FieldEntity::getEnt > > > > FieldEntity_multiIndex
MultiIndex container keeps FieldEntity.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode loop_entities(const std::string field_name, EntityMethod &method, Range const *const ents=nullptr, int verb=DEFAULT_VERBOSITY)=0
Loop over field entities.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasVector< double > VectorDouble
char FieldBitNumber
Field bit number.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr auto field_name
Data from Cubit blocksets.
FTensor::Tensor1< double, 3 > t_centr
MoFEMErrorCode getBcData(DisplacementCubitBcData &mydata, const MoFEM::CubitMeshSets *it)
MoFEMErrorCode getEntitiesFromBc(MoFEM::Interface &mField, const MoFEM::CubitMeshSets *it)
VectorDouble scaled_values
FTensor::Tensor1< double, 3 > t_normal
VectorDouble initial_values
DirichletDisplacementBc(MoFEM::Interface &m_field, const std::string &field_name, Mat Aij, Vec X, Vec F, string blockset_name="DISPLACEMENT")
const std::string blocksetName
MoFEMErrorCode postProcess()
function is run at the end of loop
std::map< DofIdx, FieldData > mapZeroRows
MoFEMErrorCode getRotationBcFromBlock(std::vector< DataFromBc > &bc_data)
Get the Rotation Bc From Block object Use ROTATION blockset name with 7 atributes: 1,...
const std::string fieldName
field name to set Dirichlet BC
double dIag
diagonal value set on zeroed column and rows
std::vector< int > dofsIndices
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode getBcDataFromSetsAndBlocks(std::vector< DataFromBc > &bc_data)
Get the Bc Data From Sets And Blocks object Use DISPLACEMENT blockset name (default) with 6 atributes...
MoFEMErrorCode applyScaleBcData(DataFromBc &bc_data)
std::vector< double > dofsXValues
boost::ptr_vector< MethodForForceScaling > methodsOp
std::vector< double > dofsValues
MoFEM::Interface & mField
virtual MoFEMErrorCode iNitialize()
MoFEMErrorCode calculateRotationForDof(VectorDouble3 &coords, DataFromBc &bc_data)
Calculate displacements from rotation for particular dof.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual boost::shared_ptr< EntityMethod > getEntMethodPtr(DataFromBc &data)
boost::shared_ptr< vector< DataFromBc > > bcDataPtr
MoFEMErrorCode iNitialize()
MoFEMErrorCode iNitialize()
std::vector< std::string > fieldNames
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode iNitialize()
std::vector< std::string > fixFields
std::string materialPositions
MoFEMErrorCode iNitialize()
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
const Problem * problemPtr
raw pointer to problem
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
this struct keeps basic methods for moab meshset about material and boundary conditions
MoFEMErrorCode getAttributes(std::vector< double > &attributes) const
get Cubit block attributes
MoFEMErrorCode getBcDataStructure(CUBIT_BC_DATA_TYPE &data) const
MoFEMErrorCode getMeshsetIdEntitiesByDimension(Interface &moab, const int dimension, Range &entities, const bool recursive=false) const
get entities form meshset
Deprecated interface functions.
Definition of the displacement bc data structure.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients.
keeps basic data about problem
auto & getNumeredColDofsPtr() const
get access to numeredColDofsPtr storing DOFs on cols
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Mat & snes_B
preconditioner of jacobian matrix
Vec & ts_F
residual vector
Mat & ts_B
Preconditioner for ts_A.
Definition of the temperature bc data structure.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEM::Interface & mField
ReactionsMap reactionsMap
MoFEMErrorCode calculateReactions(Vec &internal)
calculate reactions from a given vector