v0.14.0
SimpleRodElement.cpp
Go to the documentation of this file.
1 /* \file SimpleRodElement.cpp
2  \brief Implementation of SimpleRod element on eDges
3 */
4 
5 
6 
7 #include <MoFEM.hpp>
8 using namespace MoFEM;
9 
10 #include <SimpleRodElement.hpp>
11 using namespace boost::numeric;
12 
14  int iD;
15 
19 
21 
23  : simpleRodYoungModulus(-1), simpleRodSectionArea(-1),
24  simpleRodPreStress(-1) {}
25 };
26 
28 
29  boost::shared_ptr<MatrixDouble> gradDispPtr =
30  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
31  boost::shared_ptr<MatrixDouble> xAtPts =
32  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
33  boost::shared_ptr<MatrixDouble> xInitAtPts =
34  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
35 
39 
40  std::map<int, BlockOptionDataSimpleRods> mapSimpleRod;
41  // ~DataAtIntegrationPtsSimpleRods() {}
42  DataAtIntegrationPtsSimpleRods(MoFEM::Interface &m_field) : mField(m_field) {
43 
44  ierr = setBlocks();
45  CHKERRABORT(PETSC_COMM_WORLD, ierr);
46  }
47 
49  MoFEMFunctionBegin; // They will be overwritten by BlockData
50  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem", "none");
51 
52  ierr = PetscOptionsEnd();
53  CHKERRQ(ierr);
55  }
56 
59 
60  simpleRodYoungModulus = data.simpleRodYoungModulus;
61  simpleRodSectionArea = data.simpleRodSectionArea;
62  simpleRodPreStress = data.simpleRodPreStress;
63 
65  }
66 
69 
71  if (bit->getName().compare(0, 3, "ROD") == 0) {
72 
73  const int id = bit->getMeshsetId();
74  mapSimpleRod[id].eDges.clear();
75  CHKERR mField.get_moab().get_entities_by_type(
76  bit->getMeshset(), MBEDGE, mapSimpleRod[id].eDges, true);
77 
78  std::vector<double> attributes;
79  bit->getAttributes(attributes);
80  if (attributes.size() < 3) {
81  SETERRQ1(
82  PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
83  "Input mesh for ROD should have 3 attributes but there is %d",
84  attributes.size());
85  }
86  mapSimpleRod[id].iD = id;
87  mapSimpleRod[id].simpleRodYoungModulus = attributes[0];
88  mapSimpleRod[id].simpleRodSectionArea = attributes[1];
89  mapSimpleRod[id].simpleRodPreStress = attributes[2];
90 
91  // Print ROD blocks after being read
92  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\nSimple rod block %d\n", id);
93  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung's modulus %3.4g\n",
94  attributes[0]);
95  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tCross-section area %3.4g\n",
96  attributes[1]);
97  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPrestress %3.4g\n",
98  attributes[2]);
99  }
100  }
101 
103  }
104 
105 private:
107 };
108 
109 /** * @brief Assemble contribution of SimpleRod element to LHS *
110  *
111  */
113 
114  boost::shared_ptr<DataAtIntegrationPtsSimpleRods> commonDataPtr;
116 
119 
121  boost::shared_ptr<DataAtIntegrationPtsSimpleRods> &common_data_ptr,
122  BlockOptionDataSimpleRods &data, const std::string field_name)
124  field_name.c_str(), field_name.c_str(), OPROWCOL),
125  commonDataPtr(common_data_ptr), dAta(data) {
126  sYmm = false;
127  }
128 
129  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
130  EntityType col_type,
131  EntitiesFieldData::EntData &row_data,
132  EntitiesFieldData::EntData &col_data) {
134 
135  // check if the edge have associated degrees of freedom
136  const int row_nb_dofs = row_data.getIndices().size();
137  if (!row_nb_dofs)
139 
140  const int col_nb_dofs = col_data.getIndices().size();
141  if (!col_nb_dofs)
143 
144  if (dAta.eDges.find(getFEEntityHandle()) == dAta.eDges.end()) {
146  }
147 
148  CHKERR commonDataPtr->getBlockData(dAta);
149  // size associated to the entity
150  locK.resize(row_nb_dofs, col_nb_dofs, false);
151  locK.clear();
152 
153  double tension_stiffness = commonDataPtr->simpleRodYoungModulus *
154  commonDataPtr->simpleRodSectionArea;
155 
156  VectorDouble coords;
157  coords = getCoords();
158  double L = getLength();
159  double coeff = tension_stiffness / (L * L * L);
160 
161  double x21 = coords(3) - coords(0);
162  double y21 = coords(4) - coords(1);
163  double z21 = coords(5) - coords(2);
164 
165  // Calculate element matrix
166  locK(0, 0) = coeff * x21 * x21;
167  locK(0, 1) = coeff * x21 * y21;
168  locK(0, 2) = coeff * x21 * z21;
169  locK(0, 3) = -coeff * x21 * x21;
170  locK(0, 4) = -coeff * x21 * y21;
171  locK(0, 5) = -coeff * x21 * z21;
172 
173  locK(1, 0) = locK(0, 1);
174  locK(1, 1) = coeff * y21 * y21;
175  locK(1, 2) = coeff * y21 * z21;
176  locK(1, 3) = -coeff * y21 * x21;
177  locK(1, 4) = -coeff * y21 * y21;
178  locK(1, 5) = -coeff * y21 * z21;
179 
180  locK(2, 0) = locK(0, 2);
181  locK(2, 1) = locK(1, 2);
182  locK(2, 2) = coeff * z21 * z21;
183  locK(2, 3) = -coeff * z21 * x21;
184  locK(2, 4) = -coeff * z21 * y21;
185  locK(2, 5) = -coeff * z21 * z21;
186 
187  locK(3, 0) = locK(0, 3);
188  locK(3, 1) = locK(1, 3);
189  locK(3, 2) = locK(2, 3);
190  locK(3, 3) = coeff * x21 * x21;
191  locK(3, 4) = coeff * x21 * y21;
192  locK(3, 5) = coeff * x21 * z21;
193 
194  locK(4, 0) = locK(0, 4);
195  locK(4, 1) = locK(1, 4);
196  locK(4, 2) = locK(2, 4);
197  locK(4, 3) = locK(3, 4);
198  locK(4, 4) = coeff * y21 * y21;
199  locK(4, 5) = coeff * y21 * z21;
200 
201  locK(5, 0) = locK(0, 5);
202  locK(5, 1) = locK(1, 5);
203  locK(5, 2) = locK(2, 5);
204  locK(5, 3) = locK(3, 5);
205  locK(5, 4) = locK(4, 5);
206  locK(5, 5) = coeff * z21 * z21;
207 
208  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0), ADD_VALUES);
209 
211  }
212 };
213 
214 /** * @brief Add ROD pre-stress to the RHS *
215  */
218 
219  // vector used to store force vector for each degree of freedom
221 
222  boost::shared_ptr<DataAtIntegrationPtsSimpleRods> commonDataPtr;
224 
226  boost::shared_ptr<DataAtIntegrationPtsSimpleRods> &common_data_ptr,
227  BlockOptionDataSimpleRods &data, const std::string field_name)
229  field_name.c_str(), OPROW),
230  commonDataPtr(common_data_ptr), dAta(data) {}
231 
232  MoFEMErrorCode doWork(int side, EntityType type,
234 
236 
237  // check that the edge have associated degrees of freedom
238  const int nb_dofs = data.getIndices().size();
239  if (nb_dofs == 0)
241 
242  if (dAta.eDges.find(getFEEntityHandle()) == dAta.eDges.end()) {
244  }
245 
246  CHKERR commonDataPtr->getBlockData(dAta);
247 
248  // size of force vector associated to the entity
249  // set equal to the number of degrees of freedom of associated with the
250  // entity
251  nF.resize(nb_dofs, false);
252  nF.clear();
253 
254  double axial_force =
255  commonDataPtr->simpleRodSectionArea * commonDataPtr->simpleRodPreStress;
256 
257  auto dir = getDirection();
258  dir /= norm_2(dir);
259  for (auto d : {0, 1, 2}) {
260  nF(d) = -axial_force * dir[d];
261  nF(d + 3) = axial_force * dir[d];
262  }
263 
264  CHKERR VecSetValues(getKSPf(), data, &nF[0], ADD_VALUES);
265 
267  }
268 };
269 
271  MoFEM::Interface &m_field, const std::string field_name,
272  const std::string mesh_nodals_positions) {
274 
275  // Define boundary element that operates on rows, columns and data of a
276  // given field
277  CHKERR m_field.add_finite_element("SIMPLE_ROD", MF_ZERO);
281  if (m_field.check_field(mesh_nodals_positions)) {
282  CHKERR m_field.modify_finite_element_add_field_data("SIMPLE_ROD",
283  mesh_nodals_positions);
284  }
285  // Add entities to that element, here we add all eDges with ROD
286  // from cubit
288  if (bit->getName().compare(0, 3, "ROD") == 0) {
289  CHKERR m_field.add_ents_to_finite_element_by_type(bit->getMeshset(),
290  MBEDGE, "SIMPLE_ROD");
291  }
292  }
293  CHKERR m_field.build_finite_elements("SIMPLE_ROD");
294 
296 }
297 
299  MoFEM::Interface &m_field,
300  boost::shared_ptr<EdgeElementForcesAndSourcesCore>
301  fe_simple_rod_lhs_ptr,
302  boost::shared_ptr<EdgeElementForcesAndSourcesCore>
303  fe_simple_rod_rhs_ptr,
304  const std::string field_name, const std::string mesh_nodals_positions) {
306 
307  // Push operators to instances for SimpleRod elements
308  // loop over blocks
309  boost::shared_ptr<DataAtIntegrationPtsSimpleRods> commonDataPtr =
310  boost::make_shared<DataAtIntegrationPtsSimpleRods>(m_field);
311  CHKERR commonDataPtr->getParameters();
312 
313  for (auto &sitSimpleRod : commonDataPtr->mapSimpleRod) {
314  fe_simple_rod_lhs_ptr->getOpPtrVector().push_back(
315  new OpSimpleRodK(commonDataPtr, sitSimpleRod.second, field_name));
316 
317  fe_simple_rod_rhs_ptr->getOpPtrVector().push_back(new OpSimpleRodPreStress(
318  commonDataPtr, sitSimpleRod.second, field_name));
319  }
320 
322 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
OpSimpleRodPreStress
Add ROD pre-stress to the RHS *.
Definition: SimpleRodElement.cpp:216
DataAtIntegrationPtsSimpleRods
Definition: SimpleRodElement.cpp:27
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
BlockOptionDataSimpleRods::simpleRodSectionArea
double simpleRodSectionArea
Definition: SimpleRodElement.cpp:17
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
DataAtIntegrationPtsSimpleRods::simpleRodYoungModulus
double simpleRodYoungModulus
Definition: SimpleRodElement.cpp:36
OpSimpleRodPreStress::commonDataPtr
boost::shared_ptr< DataAtIntegrationPtsSimpleRods > commonDataPtr
Definition: SimpleRodElement.cpp:222
OpSimpleRodK::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: SimpleRodElement.cpp:129
MoFEM.hpp
OpSimpleRodK
Assemble contribution of SimpleRod element to LHS *.
Definition: SimpleRodElement.cpp:112
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
OpSimpleRodK::dAta
BlockOptionDataSimpleRods & dAta
Definition: SimpleRodElement.cpp:115
DataAtIntegrationPtsSimpleRods::getParameters
MoFEMErrorCode getParameters()
Definition: SimpleRodElement.cpp:48
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
OpSimpleRodK::commonDataPtr
boost::shared_ptr< DataAtIntegrationPtsSimpleRods > commonDataPtr
Definition: SimpleRodElement.cpp:114
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
OpSimpleRodPreStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: SimpleRodElement.cpp:232
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MetaSimpleRodElement::setSimpleRodOperators
static MoFEMErrorCode setSimpleRodOperators(MoFEM::Interface &m_field, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_lhs_ptr, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Implementation of SimpleRod element. Set operators to calculate LHS and RHS.
Definition: SimpleRodElement.cpp:298
BlockOptionDataSimpleRods::simpleRodPreStress
double simpleRodPreStress
Definition: SimpleRodElement.cpp:18
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
SimpleRodElement.hpp
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
DataAtIntegrationPtsSimpleRods::mapSimpleRod
std::map< int, BlockOptionDataSimpleRods > mapSimpleRod
Definition: SimpleRodElement.cpp:40
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
OpSimpleRodPreStress::OpSimpleRodPreStress
OpSimpleRodPreStress(boost::shared_ptr< DataAtIntegrationPtsSimpleRods > &common_data_ptr, BlockOptionDataSimpleRods &data, const std::string field_name)
Definition: SimpleRodElement.cpp:225
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
BlockOptionDataSimpleRods::iD
int iD
Definition: SimpleRodElement.cpp:14
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
DataAtIntegrationPtsSimpleRods::setBlocks
MoFEMErrorCode setBlocks()
Definition: SimpleRodElement.cpp:67
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
OpSimpleRodPreStress::dAta
BlockOptionDataSimpleRods & dAta
Definition: SimpleRodElement.cpp:223
OpSimpleRodPreStress::nF
VectorDouble nF
Definition: SimpleRodElement.cpp:220
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
DataAtIntegrationPtsSimpleRods::getBlockData
MoFEMErrorCode getBlockData(BlockOptionDataSimpleRods &data)
Definition: SimpleRodElement.cpp:57
OpSimpleRodK::OpSimpleRodK
OpSimpleRodK(boost::shared_ptr< DataAtIntegrationPtsSimpleRods > &common_data_ptr, BlockOptionDataSimpleRods &data, const std::string field_name)
Definition: SimpleRodElement.cpp:120
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
DataAtIntegrationPtsSimpleRods::simpleRodPreStress
double simpleRodPreStress
Definition: SimpleRodElement.cpp:38
DataAtIntegrationPtsSimpleRods::mField
MoFEM::Interface & mField
Definition: SimpleRodElement.cpp:106
BlockOptionDataSimpleRods::eDges
Range eDges
Definition: SimpleRodElement.cpp:20
Range
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MetaSimpleRodElement::addSimpleRodElements
static MoFEMErrorCode addSimpleRodElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare SimpleRod element.
Definition: SimpleRodElement.cpp:270
BlockOptionDataSimpleRods::simpleRodYoungModulus
double simpleRodYoungModulus
Definition: SimpleRodElement.cpp:16
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
BlockOptionDataSimpleRods
Definition: SimpleRodElement.cpp:13
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
DataAtIntegrationPtsSimpleRods::DataAtIntegrationPtsSimpleRods
DataAtIntegrationPtsSimpleRods(MoFEM::Interface &m_field)
Definition: SimpleRodElement.cpp:42
OpSimpleRodK::locK
MatrixDouble locK
Definition: SimpleRodElement.cpp:117
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
BlockOptionDataSimpleRods::BlockOptionDataSimpleRods
BlockOptionDataSimpleRods()
Definition: SimpleRodElement.cpp:22
DataAtIntegrationPtsSimpleRods::simpleRodSectionArea
double simpleRodSectionArea
Definition: SimpleRodElement.cpp:37
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
OpSimpleRodK::transLocK
MatrixDouble transLocK
Definition: SimpleRodElement.cpp:118