v0.13.1
Loading...
Searching...
No Matches
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>
8using namespace MoFEM;
9
10#include <SimpleRodElement.hpp>
11using namespace boost::numeric;
12
14 int iD;
15
19
21
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() {}
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
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
105private:
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,
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
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}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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 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 bool check_field(const std::string &name) const =0
check if field is in database
#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.
auto bit
set bit
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr auto field_name
DataAtIntegrationPtsSimpleRods(MoFEM::Interface &m_field)
MoFEMErrorCode getBlockData(BlockOptionDataSimpleRods &data)
boost::shared_ptr< MatrixDouble > xAtPts
boost::shared_ptr< MatrixDouble > gradDispPtr
boost::shared_ptr< MatrixDouble > xInitAtPts
std::map< int, BlockOptionDataSimpleRods > mapSimpleRod
static MoFEMErrorCode addSimpleRodElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare SimpleRod element.
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.
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
Assemble contribution of SimpleRod element to LHS *.
MatrixDouble transLocK
BlockOptionDataSimpleRods & dAta
boost::shared_ptr< DataAtIntegrationPtsSimpleRods > commonDataPtr
MatrixDouble locK
OpSimpleRodK(boost::shared_ptr< DataAtIntegrationPtsSimpleRods > &common_data_ptr, BlockOptionDataSimpleRods &data, const std::string field_name)
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.
Add ROD pre-stress to the RHS *.
OpSimpleRodPreStress(boost::shared_ptr< DataAtIntegrationPtsSimpleRods > &common_data_ptr, BlockOptionDataSimpleRods &data, const std::string field_name)
boost::shared_ptr< DataAtIntegrationPtsSimpleRods > commonDataPtr
BlockOptionDataSimpleRods & dAta
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.