11 using namespace boost::numeric;
23 : simpleRodYoungModulus(-1), simpleRodSectionArea(-1),
24 simpleRodPreStress(-1) {}
29 boost::shared_ptr<MatrixDouble> gradDispPtr =
31 boost::shared_ptr<MatrixDouble> xAtPts =
33 boost::shared_ptr<MatrixDouble> xInitAtPts =
45 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
50 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Problem",
"none");
52 ierr = PetscOptionsEnd();
71 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
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);
78 std::vector<double> attributes;
79 bit->getAttributes(attributes);
80 if (attributes.size() < 3) {
83 "Input mesh for ROD should have 3 attributes but there is %d",
86 mapSimpleRod[id].iD = id;
87 mapSimpleRod[id].simpleRodYoungModulus = attributes[0];
88 mapSimpleRod[id].simpleRodSectionArea = attributes[1];
89 mapSimpleRod[id].simpleRodPreStress = attributes[2];
92 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"\nSimple rod block %d\n",
id);
93 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"\tYoung's modulus %3.4g\n",
95 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"\tCross-section area %3.4g\n",
97 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"\tPrestress %3.4g\n",
121 boost::shared_ptr<DataAtIntegrationPtsSimpleRods> &common_data_ptr,
125 commonDataPtr(common_data_ptr), dAta(data) {
136 const int row_nb_dofs = row_data.
getIndices().size();
140 const int col_nb_dofs = col_data.
getIndices().size();
144 if (dAta.
eDges.find(getFEEntityHandle()) == dAta.
eDges.end()) {
148 CHKERR commonDataPtr->getBlockData(dAta);
150 locK.resize(row_nb_dofs, col_nb_dofs,
false);
153 double tension_stiffness = commonDataPtr->simpleRodYoungModulus *
154 commonDataPtr->simpleRodSectionArea;
157 coords = getCoords();
158 double L = getLength();
159 double coeff = tension_stiffness / (
L *
L *
L);
161 double x21 = coords(3) - coords(0);
162 double y21 = coords(4) - coords(1);
163 double z21 = coords(5) - coords(2);
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;
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;
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;
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;
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;
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;
226 boost::shared_ptr<DataAtIntegrationPtsSimpleRods> &common_data_ptr,
230 commonDataPtr(common_data_ptr), dAta(data) {}
242 if (dAta.
eDges.find(getFEEntityHandle()) == dAta.
eDges.end()) {
246 CHKERR commonDataPtr->getBlockData(dAta);
251 nF.resize(nb_dofs,
false);
255 commonDataPtr->simpleRodSectionArea * commonDataPtr->simpleRodPreStress;
257 auto dir = getDirection();
259 for (
auto d : {0, 1, 2}) {
260 nF(
d) = -axial_force * dir[
d];
261 nF(
d + 3) = axial_force * dir[
d];
272 const std::string mesh_nodals_positions) {
283 mesh_nodals_positions);
288 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
290 MBEDGE,
"SIMPLE_ROD");
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) {
309 boost::shared_ptr<DataAtIntegrationPtsSimpleRods> commonDataPtr =
310 boost::make_shared<DataAtIntegrationPtsSimpleRods>(m_field);
311 CHKERR commonDataPtr->getParameters();
313 for (
auto &sitSimpleRod : commonDataPtr->mapSimpleRod) {
314 fe_simple_rod_lhs_ptr->getOpPtrVector().push_back(
318 commonDataPtr, sitSimpleRod.second,
field_name));