v0.15.0
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
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 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem", "none");
51
52 PetscOptionsEnd();
54 }
55
65
68
70 if (bit->getName().compare(0, 3, "ROD") == 0) {
71
72 const int id = bit->getMeshsetId();
73 mapSimpleRod[id].eDges.clear();
74 CHKERR mField.get_moab().get_entities_by_type(
75 bit->getMeshset(), MBEDGE, mapSimpleRod[id].eDges, true);
76
77 std::vector<double> attributes;
78 bit->getAttributes(attributes);
79 if (attributes.size() < 3) {
80 SETERRQ(
81 PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
82 "Input mesh for ROD should have 3 attributes but there is %ld",
83 attributes.size());
84 }
85 mapSimpleRod[id].iD = id;
86 mapSimpleRod[id].simpleRodYoungModulus = attributes[0];
87 mapSimpleRod[id].simpleRodSectionArea = attributes[1];
88 mapSimpleRod[id].simpleRodPreStress = attributes[2];
89
90 // Print ROD blocks after being read
91 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\nSimple rod block %d\n", id);
92 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung's modulus %3.4g\n",
93 attributes[0]);
94 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tCross-section area %3.4g\n",
95 attributes[1]);
96 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPrestress %3.4g\n",
97 attributes[2]);
98 }
99 }
100
102 }
103
104private:
106};
107
108/** * @brief Assemble contribution of SimpleRod element to LHS *
109 *
110 */
112
113 boost::shared_ptr<DataAtIntegrationPtsSimpleRods> commonDataPtr;
115
118
120 boost::shared_ptr<DataAtIntegrationPtsSimpleRods> &common_data_ptr,
121 BlockOptionDataSimpleRods &data, const std::string field_name)
123 field_name.c_str(), field_name.c_str(), OPROWCOL),
124 commonDataPtr(common_data_ptr), dAta(data) {
125 sYmm = false;
126 }
127
128 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
129 EntityType col_type,
131 EntitiesFieldData::EntData &col_data) {
133
134 // check if the edge have associated degrees of freedom
135 const int row_nb_dofs = row_data.getIndices().size();
136 if (!row_nb_dofs)
138
139 const int col_nb_dofs = col_data.getIndices().size();
140 if (!col_nb_dofs)
142
143 if (dAta.eDges.find(getFEEntityHandle()) == dAta.eDges.end()) {
145 }
146
147 CHKERR commonDataPtr->getBlockData(dAta);
148 // size associated to the entity
149 locK.resize(row_nb_dofs, col_nb_dofs, false);
150 locK.clear();
151
152 double tension_stiffness = commonDataPtr->simpleRodYoungModulus *
153 commonDataPtr->simpleRodSectionArea;
154
155 VectorDouble coords;
156 coords = getCoords();
157 double L = getLength();
158 double coeff = tension_stiffness / (L * L * L);
159
160 double x21 = coords(3) - coords(0);
161 double y21 = coords(4) - coords(1);
162 double z21 = coords(5) - coords(2);
163
164 // Calculate element matrix
165 locK(0, 0) = coeff * x21 * x21;
166 locK(0, 1) = coeff * x21 * y21;
167 locK(0, 2) = coeff * x21 * z21;
168 locK(0, 3) = -coeff * x21 * x21;
169 locK(0, 4) = -coeff * x21 * y21;
170 locK(0, 5) = -coeff * x21 * z21;
171
172 locK(1, 0) = locK(0, 1);
173 locK(1, 1) = coeff * y21 * y21;
174 locK(1, 2) = coeff * y21 * z21;
175 locK(1, 3) = -coeff * y21 * x21;
176 locK(1, 4) = -coeff * y21 * y21;
177 locK(1, 5) = -coeff * y21 * z21;
178
179 locK(2, 0) = locK(0, 2);
180 locK(2, 1) = locK(1, 2);
181 locK(2, 2) = coeff * z21 * z21;
182 locK(2, 3) = -coeff * z21 * x21;
183 locK(2, 4) = -coeff * z21 * y21;
184 locK(2, 5) = -coeff * z21 * z21;
185
186 locK(3, 0) = locK(0, 3);
187 locK(3, 1) = locK(1, 3);
188 locK(3, 2) = locK(2, 3);
189 locK(3, 3) = coeff * x21 * x21;
190 locK(3, 4) = coeff * x21 * y21;
191 locK(3, 5) = coeff * x21 * z21;
192
193 locK(4, 0) = locK(0, 4);
194 locK(4, 1) = locK(1, 4);
195 locK(4, 2) = locK(2, 4);
196 locK(4, 3) = locK(3, 4);
197 locK(4, 4) = coeff * y21 * y21;
198 locK(4, 5) = coeff * y21 * z21;
199
200 locK(5, 0) = locK(0, 5);
201 locK(5, 1) = locK(1, 5);
202 locK(5, 2) = locK(2, 5);
203 locK(5, 3) = locK(3, 5);
204 locK(5, 4) = locK(4, 5);
205 locK(5, 5) = coeff * z21 * z21;
206
207 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0), ADD_VALUES);
208
210 }
211};
212
213/** * @brief Add ROD pre-stress to the RHS *
214 */
217
218 // vector used to store force vector for each degree of freedom
220
221 boost::shared_ptr<DataAtIntegrationPtsSimpleRods> commonDataPtr;
223
225 boost::shared_ptr<DataAtIntegrationPtsSimpleRods> &common_data_ptr,
226 BlockOptionDataSimpleRods &data, const std::string field_name)
228 field_name.c_str(), OPROW),
229 commonDataPtr(common_data_ptr), dAta(data) {}
230
231 MoFEMErrorCode doWork(int side, EntityType type,
233
235
236 // check that the edge have associated degrees of freedom
237 const int nb_dofs = data.getIndices().size();
238 if (nb_dofs == 0)
240
241 if (dAta.eDges.find(getFEEntityHandle()) == dAta.eDges.end()) {
243 }
244
245 CHKERR commonDataPtr->getBlockData(dAta);
246
247 // size of force vector associated to the entity
248 // set equal to the number of degrees of freedom of associated with the
249 // entity
250 nF.resize(nb_dofs, false);
251 nF.clear();
252
253 double axial_force =
254 commonDataPtr->simpleRodSectionArea * commonDataPtr->simpleRodPreStress;
255
256 auto dir = getDirection();
257 dir /= norm_2(dir);
258 for (auto d : {0, 1, 2}) {
259 nF(d) = -axial_force * dir[d];
260 nF(d + 3) = axial_force * dir[d];
261 }
262
263 CHKERR VecSetValues(getKSPf(), data, &nF[0], ADD_VALUES);
264
266 }
267};
268
270 MoFEM::Interface &m_field, const std::string field_name,
271 const std::string mesh_nodals_positions) {
273
274 // Define boundary element that operates on rows, columns and data of a
275 // given field
276 CHKERR m_field.add_finite_element("SIMPLE_ROD", MF_ZERO);
280 if (m_field.check_field(mesh_nodals_positions)) {
281 CHKERR m_field.modify_finite_element_add_field_data("SIMPLE_ROD",
282 mesh_nodals_positions);
283 }
284 // Add entities to that element, here we add all eDges with ROD
285 // from cubit
287 if (bit->getName().compare(0, 3, "ROD") == 0) {
288 CHKERR m_field.add_ents_to_finite_element_by_type(bit->getMeshset(),
289 MBEDGE, "SIMPLE_ROD");
290 }
291 }
292 CHKERR m_field.build_finite_elements("SIMPLE_ROD");
293
295}
296
298 MoFEM::Interface &m_field,
299 boost::shared_ptr<EdgeElementForcesAndSourcesCore>
300 fe_simple_rod_lhs_ptr,
301 boost::shared_ptr<EdgeElementForcesAndSourcesCore>
302 fe_simple_rod_rhs_ptr,
303 const std::string field_name, const std::string mesh_nodals_positions) {
305
306 // Push operators to instances for SimpleRod elements
307 // loop over blocks
308 boost::shared_ptr<DataAtIntegrationPtsSimpleRods> commonDataPtr =
309 boost::make_shared<DataAtIntegrationPtsSimpleRods>(m_field);
310 CHKERR commonDataPtr->getParameters();
311
312 for (auto &sitSimpleRod : commonDataPtr->mapSimpleRod) {
313 fe_simple_rod_lhs_ptr->getOpPtrVector().push_back(
314 new OpSimpleRodK(commonDataPtr, sitSimpleRod.second, field_name));
315
316 fe_simple_rod_rhs_ptr->getOpPtrVector().push_back(new OpSimpleRodPreStress(
317 commonDataPtr, sitSimpleRod.second, field_name));
318 }
319
321}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
@ MF_ZERO
#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 ...
@ BLOCKSET
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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 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
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 modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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.