v0.10.0
ThermalStressElement.hpp
Go to the documentation of this file.
1 /** \file ThermalStressElement.hpp
2  \ingroup mofem_thermal_elem
3  \brief Implementation of thermal stresses element
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #ifndef __THERMALSTRESSELEMENT_HPP
22 #define __THERMALSTRESSELEMENT_HPP
23 
24 /** \brief Implentation of thermal stress element
25 \ingroup mofem_thermal_elem
26 */
28 
32  int getRule(int order) { return 2 * (order - 1); };
33  };
34 
37 
40  : feThermalStressRhs(m_field), mField(m_field) {}
41 
42  struct BlockData {
43  double youngModulus;
44  double poissonRatio;
48  Range tEts;
49  };
50  std::map<int, BlockData> setOfBlocks;
51 
52  struct CommonData {
54  };
56 
59 
61  int verb;
62  OpGetTemperatureAtGaussPts(const std::string field_name,
63  CommonData &common_data, int _verb = 0)
65  field_name, UserDataOperator::OPROW),
66  commonData(common_data), verb(_verb) {}
67 
68  MoFEMErrorCode doWork(int side, EntityType type,
71  if (data.getFieldData().size() == 0)
73  int nb_dofs = data.getFieldData().size();
74  int nb_gauss_pts = data.getN().size1();
75  // initialize
76  commonData.temperatureAtGaussPts.resize(nb_gauss_pts);
77  if (type == MBVERTEX) {
78  std::fill(commonData.temperatureAtGaussPts.begin(),
80  }
81  for (int gg = 0; gg < nb_gauss_pts; gg++) {
83  inner_prod(data.getN(gg, nb_dofs), data.getFieldData());
84  }
86  }
87  };
88 
91 
92  Vec F;
95  int verb;
96  OpThermalStressRhs(const std::string field_name, Vec _F, BlockData &data,
97  CommonData &common_data, int _verb = 0)
99  field_name, UserDataOperator::OPROW),
100  F(_F), dAta(data), commonData(common_data), verb(_verb) {}
101 
103  MoFEMErrorCode doWork(int side, EntityType type,
106 
107  if (data.getIndices().size() == 0)
109  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
110  dAta.tEts.end())
112 
113  const auto &dof_ptr = data.getFieldDofs()[0];
114  int rank = dof_ptr->getNbOfCoeffs();
115  if (rank != 3) {
116  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
117  "data inconsistency");
118  }
119 
120  unsigned int nb_dofs = data.getIndices().size();
121  if (nb_dofs % rank != 0) {
122  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
123  "data inconsistency");
124  }
125  if (data.getN().size2() < nb_dofs / rank) {
126  SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
127  "data inconsistency N.size2 %d nb_dof %d rank %d",
128  data.getN().size2(), nb_dofs, rank);
129  }
130 
131  Nf.resize(nb_dofs);
132  bzero(&*Nf.data().begin(), nb_dofs * sizeof(FieldData));
133 
134  if (verb > VERBOSE) {
135  if (type == MBVERTEX) {
136  cout << commonData.temperatureAtGaussPts << std::endl;
137  cout << "thermal expansion " << dAta.thermalExpansion << std::endl;
138  }
139  }
140 
141  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
142 
144  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "invalid data");
145  }
146 
147  double phi =
149  double eps_thermal = dAta.thermalExpansion * phi;
150  double sig_thermal =
151  -eps_thermal * (dAta.youngModulus / (1. - 2 * dAta.poissonRatio));
152 
153  double val = sig_thermal * getVolume() * getGaussPts()(3, gg);
154 
155  double *diff_N;
156  diff_N = &data.getDiffN()(gg, 0);
157  cblas_daxpy(nb_dofs, val, diff_N, 1, &Nf[0], 1);
158  }
159 
160  CHKERR VecSetValues(F, data.getIndices().size(), &data.getIndices()[0],
161  &Nf[0], ADD_VALUES);
162 
164  }
165  };
166 
168  const std::string fe_name, const std::string field_name,
169  const std::string thermal_field_name,
170  const std::string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
172  if (mField.check_field(thermal_field_name)) {
173 
179  thermal_field_name);
180  if (mField.check_field(mesh_nodals_positions)) {
182  fe_name, mesh_nodals_positions);
183  }
185  mField, BLOCKSET | MAT_ELASTICSET, it)) {
186  Mat_Elastic mydata;
187  CHKERR it->getAttributeDataStructure(mydata);
188  setOfBlocks[it->getMeshsetId()].youngModulus = mydata.data.Young;
189  setOfBlocks[it->getMeshsetId()].poissonRatio = mydata.data.Poisson;
190  setOfBlocks[it->getMeshsetId()].thermalExpansion =
191  mydata.data.ThermalExpansion;
192  CHKERR mField.get_moab().get_entities_by_type(
193  it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
195  setOfBlocks[it->getMeshsetId()].tEts, MBTET, fe_name);
196  double ref_temp;
197  PetscBool flg;
198  CHKERR PetscOptionsGetReal(PETSC_NULL, PETSC_NULL, "-my_ref_temp",
199  &ref_temp, &flg);
200  if (flg == PETSC_TRUE) {
201  PetscPrintf(mField.get_comm(), "set refernce temperature %3.2f\n",
202  ref_temp);
203  setOfBlocks[it->getMeshsetId()].refTemperature = ref_temp;
204  }
205  }
206  }
208  }
209 
210  /// \deprecated Do not use this fiction with spelling mistake
212  const std::string fe_name, const std::string field_name,
213  const std::string thermal_field_name,
214  const std::string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
215  return addThermalStressElement(fe_name, field_name, thermal_field_name,
216  mesh_nodals_positions);
217  }
218 
220  string thermal_field_name, Vec &F,
221  int verb = 0) {
223  if (mField.check_field(thermal_field_name)) {
224  std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
225  for (; sit != setOfBlocks.end(); sit++) {
226  // add finite elemen
228  new OpGetTemperatureAtGaussPts(thermal_field_name, commonData,
229  verb));
231  field_name, F, sit->second, commonData, verb));
232  }
233  }
235  }
236 };
237 
238 #endif //__THERMALSTRESSELEMENT_HPP
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
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 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
Deprecated interface functions.
Implentation of thermal stress element.
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
virtual moab::Interface & get_moab()=0
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
OpGetTemperatureAtGaussPts(const std::string field_name, CommonData &common_data, int _verb=0)
MoFEMErrorCode addThermalStressElement(const std::string fe_name, const std::string field_name, const std::string thermal_field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:604
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
block name is "MAT_ELASTIC"
Definition: definitions.h:228
double FieldData
Field data type.
Definition: Types.hpp:36
constexpr int order
#define DEPRECATED
Definition: definitions.h:29
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
virtual bool check_field(const std::string &name) const =0
check if field is in database
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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
MoFEM::Interface & mField
virtual MPI_Comm & get_comm() const =0
MoFEMErrorCode setThermalStressRhsOperators(string field_name, string thermal_field_name, Vec &F, int verb=0)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
#define _F(n)
Definition: quad.c:25
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
ThermalStressElement(MoFEM::Interface &m_field)
OpThermalStressRhs(const std::string field_name, Vec _F, BlockData &data, CommonData &common_data, int _verb=0)
MyVolumeFE(MoFEM::Interface &m_field)
std::map< int, BlockData > setOfBlocks
DEPRECATED MoFEMErrorCode addThermalSterssElement(const std::string fe_name, const std::string field_name, const std::string thermal_field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
MyVolumeFE & getLoopThermalStressRhs()