1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18 // <h1>Miscellaneous Example 12 - MITC4 Shell Elements</h1>
19 // \author Sylvain Vallaghe
20 // \date 2016
21 //
22 // This example implements a variation of MITC4 shell elements.
23 // The infamous "pinched cylinder" problem is solved and the solution
24 // is compared to analytical values at selected points.
25 // The implementation follows very closely the notations of the french
26 // reference book on structural analysis with finite elements:
27 // Batoz, Jean-Louis, and Gouri Dhatt.
28 // "Modelisation des structures par elements finis, Vol. 3: Coques." Hermes, Paris (1992).
29
30 // C++ include files that we need
31 #include <iostream>
32
33 // LibMesh includes
34 #include "libmesh/libmesh.h"
35 #include "libmesh/mesh.h"
36 #include "libmesh/linear_implicit_system.h"
37 #include "libmesh/equation_systems.h"
38 #include "libmesh/fe.h"
39 #include "libmesh/quadrature.h"
40 #include "libmesh/node.h"
41 #include "libmesh/elem.h"
42 #include "libmesh/dof_map.h"
43 #include "libmesh/quadrature_gauss.h"
44 #include "libmesh/vector_value.h"
45 #include "libmesh/tensor_value.h"
46 #include "libmesh/dense_matrix.h"
47 #include "libmesh/dense_submatrix.h"
48 #include "libmesh/dense_vector.h"
49 #include "libmesh/dense_subvector.h"
50 #include "libmesh/sparse_matrix.h"
51 #include "libmesh/numeric_vector.h"
52 #include "libmesh/exodusII_io.h"
53 #include "libmesh/dirichlet_boundaries.h"
54 #include "libmesh/zero_function.h"
55 #include "libmesh/linear_solver.h"
56 #include "libmesh/getpot.h"
57 #include "libmesh/enum_solver_package.h"
58 #include "libmesh/enum_solver_type.h"
59 #include "libmesh/parallel.h"
60
61 // Eigen includes
62 #ifdef LIBMESH_HAVE_EIGEN
63 #include "libmesh/ignore_warnings.h"
64 # include <Eigen/Dense>
65 #include "libmesh/restore_warnings.h"
66
67 // Use Eigen dense numerics with libMesh::Real
68 typedef Eigen::Matrix<libMesh::Real, 2, 2> MyMatrix2d;
69 typedef Eigen::Matrix<libMesh::Real, 3, 3> MyMatrix3d;
70 typedef Eigen::Matrix<libMesh::Real, 2, 1> MyVector2d;
71 typedef Eigen::Matrix<libMesh::Real, 3, 1> MyVector3d;
72 typedef Eigen::Matrix<libMesh::Real, Eigen::Dynamic, Eigen::Dynamic> MyMatrixXd;
73 #endif
74
75 // Bring in everything from the libMesh namespace
76 using namespace libMesh;
77
78 // Function prototype. This is the function that will assemble
79 // the stiffness matrix and the right-hand-side vector ready
80 // for solution.
81 void assemble_shell (EquationSystems & es,
82 const std::string & system_name);
83
84 // Begin the main program.
main(int argc,char ** argv)85 int main (int argc, char ** argv)
86 {
87 // Initialize libMesh.
88 LibMeshInit init (argc, argv);
89
90 // Skip this 3D example if libMesh was compiled as 1D/2D-only.
91 libmesh_example_requires (3 == LIBMESH_DIM, "3D support");
92
93 // We use Dirichlet boundary conditions here
94 #ifndef LIBMESH_ENABLE_DIRICHLET
95 libmesh_example_requires(false, "--enable-dirichlet");
96 #endif
97
98 // Our input mesh here is in ExodusII format
99 #ifndef LIBMESH_HAVE_EXODUS_API
100 libmesh_example_requires (false, "ExodusII support");
101 #endif
102
103 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
104 libmesh_example_requires (false, "second derivatives enabled");
105 #endif
106
107 // This example does a bunch of linear algebra during assembly, and
108 // therefore requires Eigen.
109 #ifndef LIBMESH_HAVE_EIGEN
110 libmesh_example_requires(false, "--enable-eigen");
111 #endif
112
113 // This example converts between ExodusII and XDR files, therefore
114 // it requires XDR support in libmesh.
115 #ifndef LIBMESH_HAVE_XDR
116 libmesh_example_requires (false, "XDR support");
117 #endif
118
119 // This examples requires parallel minloc() support, which we don't
120 // implement yet for float128
121 #ifdef LIBMESH_DEFAULT_QUADRUPLE_PRECISION
122 libmesh_example_requires (false, "--disable-quadruple-precision");
123 #endif
124
125 // Read the "distributed_load" flag from the command line
126 GetPot command_line (argc, argv);
127 int distributed_load = 0;
128 if (command_line.search(1, "-distributed_load"))
129 distributed_load = command_line.next(distributed_load);
130
131 {
132 Mesh mesh (init.comm(), 3);
133
134 // To confirm that both ExodusII and Xdr formats work for shell
135 // meshes, we read in cylinder.exo, then write out cylinder.xdr,
136 // then read in cylinder.exo again below and use that for the rest
137 // of the example.
138 mesh.read("cylinder.exo");
139 mesh.write("cylinder.xdr");
140 }
141 Mesh mesh (init.comm(), 3);
142 mesh.read("cylinder.xdr");
143
144 // Print information about the mesh to the screen.
145 mesh.print_info();
146
147 // Create an equation systems object.
148 EquationSystems equation_systems (mesh);
149
150 // Declare the system and its variables.
151 // Create a linear implicit system named "Shell".
152 LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem> ("Shell");
153
154 // Add the three displacement variables "u", "v", "w",
155 // and the three rotational variables "theta_x", "theta_y", "theta_z".
156 // All variables are Q1 (first order on a quad mesh).
157 system.add_variable ("u");
158 system.add_variable ("v");
159 system.add_variable ("w");
160 system.add_variable ("theta_x");
161 system.add_variable ("theta_y");
162 system.add_variable ("theta_z");
163
164 // Give the system a pointer to the matrix and rhs assembly
165 // function.
166 system.attach_assemble_function (assemble_shell);
167
168 // Use the parameters of the equation systems object to
169 // tell the shell system about the material properties, the
170 // shell thickness, and the external load.
171 const Real h = 0.03;
172 const Real E = 3e10;
173 const Real nu = 0.3;
174 const Real q = 1;
175 equation_systems.parameters.set<Real> ("thickness") = h;
176 equation_systems.parameters.set<Real> ("young's modulus") = E;
177 equation_systems.parameters.set<Real> ("poisson ratio") = nu;
178 equation_systems.parameters.set<Real> ("point load") = q;
179 equation_systems.parameters.set<bool>("distributed load") = (distributed_load != 0);
180
181 // Dirichlet conditions for the pinched cylinder problem.
182 // Only one 8th of the cylinder is considered using symmetry considerations.
183 // The cylinder longitudinal axis is the y-axis.
184 // The four corners of the surface are named A(3,0,0), B(3,3,0), C(0,3,3), D(0,0,3).
185 // The point load (pinch) is applied at C in the -z direction.
186 // Edge AD is the actual edge of the cylinder and is rigid in the xz-plane.
187 // Other edges have symmetric boundary conditions.
188
189 #ifdef LIBMESH_ENABLE_DIRICHLET
190 // AB w, theta_x, theta_y
191 {
192 std::set<boundary_id_type> boundary_ids;
193 boundary_ids.insert(7);
194 unsigned int variables[] = {2, 3, 4};
195 ZeroFunction<> zf;
196
197 // Most DirichletBoundary users will want to supply a "locally
198 // indexed" functor
199 DirichletBoundary dirichlet_bc
200 (boundary_ids,
201 std::vector<unsigned int>(variables, variables+3), zf,
202 LOCAL_VARIABLE_ORDER);
203 system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
204 }
205 // BC v, theta_x, theta_z
206 {
207 std::set<boundary_id_type> boundary_ids;
208 boundary_ids.insert(8);
209 unsigned int variables[] = {1, 3, 5};
210 ZeroFunction<> zf;
211
212 DirichletBoundary dirichlet_bc
213 (boundary_ids,
214 std::vector<unsigned int>(variables, variables+3), zf,
215 LOCAL_VARIABLE_ORDER);
216 system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
217 }
218 // CD u, theta_y, theta_z
219 {
220 std::set<boundary_id_type> boundary_ids;
221 boundary_ids.insert(9);
222 unsigned int variables[] = {0, 4, 5};
223 ZeroFunction<> zf;
224
225 DirichletBoundary dirichlet_bc
226 (boundary_ids,
227 std::vector<unsigned int>(variables, variables+3), zf,
228 LOCAL_VARIABLE_ORDER);
229 system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
230 }
231 // AD u, w, theta_y
232 {
233 std::set<boundary_id_type> boundary_ids;
234 boundary_ids.insert(10);
235 unsigned int variables[] = {0, 2, 4};
236 ZeroFunction<> zf;
237
238 DirichletBoundary dirichlet_bc
239 (boundary_ids,
240 std::vector<unsigned int>(variables, variables+3), zf,
241 LOCAL_VARIABLE_ORDER);
242 system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
243 }
244 #endif // LIBMESH_ENABLE_DIRICHLET
245
246 // Initialize the data structures for the equation system.
247 equation_systems.init();
248
249 // Print information about the system to the screen.
250 equation_systems.print_info();
251
252 // This example can be run with EigenSparseLinearSolvers, but it
253 // only works with either the CG or SPARSELU types, and SparseLU
254 // turns out to be faster.
255 if (libMesh::default_solver_package() == EIGEN_SOLVERS)
256 system.get_linear_solver()->set_solver_type(SPARSELU);
257
258 // Solve the linear system.
259 system.solve();
260
261 // After solving the system, write the solution to an
262 // ExodusII output file ready for import in, e.g.,
263 // Paraview.
264 ExodusII_IO(mesh).write_equation_systems ("out.e", equation_systems);
265
266 // Compare with analytical solution for point load
267 if (distributed_load==0)
268 {
269 // Find the node nearest point C.
270 Node * node_C = nullptr;
271 Point point_C(0, 3, 3);
272 {
273 Real nearest_dist_sq = std::numeric_limits<Real>::max();
274
275 // Find the closest local node. On a DistributedMesh we may
276 // not even know about the existence of closer non-local
277 // nodes.
278 for (auto & node : mesh.local_node_ptr_range())
279 {
280 const Real dist_sq = (*node - point_C).norm_sq();
281 if (dist_sq < nearest_dist_sq)
282 {
283 nearest_dist_sq = dist_sq;
284 node_C = node;
285 }
286 }
287
288 // Check with other processors to see if any found a closer node
289 unsigned int minrank = 0;
290 system.comm().minloc(nearest_dist_sq, minrank);
291
292 // Broadcast the ID of the closest node, so every processor can
293 // see for certain whether they have it or not.
294 dof_id_type nearest_node_id = 0;
295 if (system.processor_id() == minrank)
296 nearest_node_id = node_C->id();
297 system.comm().broadcast(nearest_node_id, minrank);
298 node_C = mesh.query_node_ptr(nearest_node_id);
299 }
300
301 // Evaluate the z-displacement "w" at the node nearest C.
302 Number w = 0;
303
304 // If we know about the closest node, and if we also own the DoFs
305 // on that node, then we can evaluate the solution at that node.
306 if (node_C)
307 {
308 const unsigned int w_var = system.variable_number ("w");
309 dof_id_type w_dof = node_C->dof_number (system.number(), w_var, 0);
310 if (w_dof >= system.get_dof_map().first_dof() &&
311 w_dof < system.get_dof_map().end_dof())
312 w = system.current_solution(w_dof);
313 }
314 system.comm().sum(w);
315
316
317 Number w_C_bar = -E*h*w/q;
318 const Real w_C_bar_analytic = 164.24;
319
320 // Print the finite element solution and the analytic
321 // prediction to the screen.
322 libMesh::out << "z-displacement of the point C: " << w_C_bar << std::endl;
323 libMesh::out << "Analytic solution: " << w_C_bar_analytic << std::endl;
324
325 // Evaluate the y-displacement "v" at point D. This time we'll
326 // evaluate at the exact point, not just the closest node.
327 Point point_D(0, 0, 3);
328 const unsigned int v_var = system.variable_number ("v");
329 Number v = system.point_value(v_var, point_D);
330
331 Number v_D_bar = E*h*v/q;
332 const Real v_D_bar_analytic = 4.114;
333
334 // Print the finite element solution and the analytic
335 // prediction to the screen.
336 libMesh::out << "y-displacement of the point D: " << v_D_bar << std::endl;
337 libMesh::out << "Analytic solution: " << v_D_bar_analytic << std::endl;
338 }
339
340 // All done.
341 return 0;
342 }
343
344
345
346 // We now define the matrix and rhs vector assembly function
347 // for the shell system.
assemble_shell(EquationSystems & es,const std::string & system_name)348 void assemble_shell (EquationSystems & es,
349 const std::string & system_name)
350 {
351 // This example requires Eigen to actually work, but we should still
352 // let it compile and throw a runtime error if you don't.
353
354 // The same holds for second derivatives,
355 // since they are class-members only depending on the config.
356 #if defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
357 // It is a good idea to make sure we are assembling
358 // the proper system.
359 libmesh_assert_equal_to (system_name, "Shell");
360
361 // Get a constant reference to the mesh object.
362 const MeshBase & mesh = es.get_mesh();
363 const unsigned int dim = mesh.mesh_dimension();
364
365 // Get a reference to the shell system object.
366 LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> (system_name);
367
368 // Get the shell parameters that we need during assembly.
369 const Real h = es.parameters.get<Real> ("thickness");
370 const Real E = es.parameters.get<Real> ("young's modulus");
371 const Real nu = es.parameters.get<Real> ("poisson ratio");
372 const Real q = es.parameters.get<Real> ("point load");
373 const bool distributed_load = es.parameters.get<bool> ("distributed load");
374
375 // The membrane elastic matrix.
376 MyMatrix3d Hm;
377 Hm <<
378 1., nu, 0.,
379 nu, 1., 0.,
380 0., 0., 0.5 * (1-nu);
381 Hm *= h * E/(1-nu*nu);
382
383 // The bending elastic matrix.
384 MyMatrix3d Hf;
385 Hf <<
386 1., nu, 0.,
387 nu, 1., 0.,
388 0., 0., 0.5 * (1-nu);
389 Hf *= h*h*h/12 * E/(1-nu*nu);
390
391 // The shear elastic matrices.
392 MyMatrix2d Hc0 = MyMatrix2d::Identity();
393 Hc0 *= h * 5./6*E/(2*(1+nu));
394
395 MyMatrix2d Hc1 = MyMatrix2d::Identity();
396 Hc1 *= h*h*h/12 * 5./6*E/(2*(1+nu));
397
398 // Get the Finite Element type, this will be
399 // the same for all variables.
400 FEType fe_type = system.variable_type (0);
401
402 std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
403 QGauss qrule (dim, fe_type.default_quadrature_order());
404 fe->attach_quadrature_rule (&qrule);
405
406 // The element Jacobian * quadrature weight at each integration point.
407 const std::vector<Real> & JxW = fe->get_JxW();
408
409 // The element shape function and its derivatives evaluated at the
410 // quadrature points.
411 const std::vector<RealGradient> & dxyzdxi = fe->get_dxyzdxi();
412 const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();
413
414 const std::vector<RealGradient> & d2xyzdxi2 = fe->get_d2xyzdxi2();
415 const std::vector<RealGradient> & d2xyzdeta2 = fe->get_d2xyzdeta2();
416 const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();
417 const std::vector<std::vector<Real>> & dphidxi = fe->get_dphidxi();
418 const std::vector<std::vector<Real>> & dphideta = fe->get_dphideta();
419 const std::vector<std::vector<Real>> & phi = fe->get_phi();
420
421 // A reference to the DofMap object for this system. The DofMap
422 // object handles the index translation from node and element numbers
423 // to degree of freedom numbers.
424 const DofMap & dof_map = system.get_dof_map();
425
426 // Define data structures to contain the element stiffness matrix.
427 DenseMatrix<Number> Ke;
428 DenseSubMatrix<Number> Ke_var[6][6] =
429 {
430 {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
431 DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)},
432 {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
433 DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)},
434 {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
435 DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)},
436 {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
437 DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)},
438 {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
439 DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)},
440 {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
441 DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)}
442 };
443
444 // Define data structures to contain the element rhs vector.
445 DenseVector<Number> Fe;
446 DenseSubVector<Number> Fe_w(Fe);
447
448 std::vector<dof_id_type> dof_indices;
449 std::vector<std::vector<dof_id_type>> dof_indices_var(6);
450
451 // Now we will loop over all the elements in the mesh. We will
452 // compute the element matrix and right-hand-side contribution.
453 for (const auto & elem : mesh.active_local_element_ptr_range())
454 {
455 dof_map.dof_indices (elem, dof_indices);
456 for (unsigned int var=0; var<6; var++)
457 dof_map.dof_indices (elem, dof_indices_var[var], var);
458
459 const unsigned int n_dofs = dof_indices.size();
460 const unsigned int n_var_dofs = dof_indices_var[0].size();
461
462 // First compute element data at the nodes
463 std::vector<Point> nodes;
464 for (auto i : elem->node_index_range())
465 nodes.push_back(elem->reference_elem()->node_ref(i));
466 fe->reinit (elem, &nodes);
467
468 // Convenient notation for the element node positions
469 MyVector3d X1(elem->node_ref(0)(0), elem->node_ref(0)(1), elem->node_ref(0)(2));
470 MyVector3d X2(elem->node_ref(1)(0), elem->node_ref(1)(1), elem->node_ref(1)(2));
471 MyVector3d X3(elem->node_ref(2)(0), elem->node_ref(2)(1), elem->node_ref(2)(2));
472 MyVector3d X4(elem->node_ref(3)(0), elem->node_ref(3)(1), elem->node_ref(3)(2));
473
474 //Store covariant basis and local orthonormal basis at the nodes
475 std::vector<MyMatrix3d> F0node;
476 std::vector<MyMatrix3d> Qnode;
477 for (auto i : elem->node_index_range())
478 {
479 MyVector3d a1;
480 a1 << dxyzdxi[i](0), dxyzdxi[i](1), dxyzdxi[i](2);
481 MyVector3d a2;
482 a2 << dxyzdeta[i](0), dxyzdeta[i](1), dxyzdeta[i](2);
483 MyVector3d n;
484 n = a1.cross(a2);
485 n /= n.norm();
486 MyMatrix3d F0;
487 F0 <<
488 a1(0), a2(0), n(0),
489 a1(1), a2(1), n(1),
490 a1(2), a2(2), n(2);
491 F0node.push_back(F0);
492
493 Real nx = n(0);
494 Real ny = n(1);
495 Real C = n(2);
496 if (std::abs(1.+C)<1e-6)
497 {
498 MyMatrix3d Q;
499 Q <<
500 1, 0, 0,
501 0, -1, 0,
502 0, 0, -1;
503 Qnode.push_back(Q);
504 }
505 else
506 {
507 MyMatrix3d Q;
508 Q <<
509 C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
510 -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
511 -nx, -ny, C;
512 Qnode.push_back(Q);
513 }
514 }
515
516 Ke.resize (n_dofs, n_dofs);
517 for (unsigned int var_i=0; var_i<6; var_i++)
518 for (unsigned int var_j=0; var_j<6; var_j++)
519 Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
520
521 Fe.resize(n_dofs);
522 Fe_w.reposition(2*n_var_dofs,n_var_dofs);
523
524 // Reinit element data at the regular Gauss quadrature points
525 fe->reinit (elem);
526
527 // Now we will build the element matrix and right-hand-side.
528 for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
529 {
530
531 //Covariant basis at the quadrature point
532 MyVector3d a1;
533 a1 << dxyzdxi[qp](0), dxyzdxi[qp](1), dxyzdxi[qp](2);
534 MyVector3d a2;
535 a2 << dxyzdeta[qp](0), dxyzdeta[qp](1), dxyzdeta[qp](2);
536 MyVector3d n;
537 n = a1.cross(a2);
538 n /= n.norm();
539 MyMatrix3d F0;
540 F0 <<
541 a1(0), a2(0), n(0),
542 a1(1), a2(1), n(1),
543 a1(2), a2(2), n(2);
544
545 //Contravariant basis
546 MyMatrix3d F0it;
547 F0it = F0.inverse().transpose();
548
549 //Local orthonormal basis at the quadrature point
550 Real nx = n(0);
551 Real ny = n(1);
552 Real C = n(2);
553 MyMatrix3d Q;
554 if (std::abs(1.+C) < 1e-6)
555 {
556 Q <<
557 1, 0, 0,
558 0, -1, 0,
559 0, 0, -1;
560 }
561 else
562 {
563 Q <<
564 C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
565 -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
566 -nx, -ny, C;
567 }
568
569 MyMatrix2d C0;
570 C0 = F0it.block<3,2>(0,0).transpose()*Q.block<3,2>(0,0);
571
572 // Normal derivatives in reference coordinates
573 MyVector3d d2Xdxi2(d2xyzdxi2[qp](0), d2xyzdxi2[qp](1), d2xyzdxi2[qp](2));
574 MyVector3d d2Xdeta2(d2xyzdeta2[qp](0), d2xyzdeta2[qp](1), d2xyzdeta2[qp](2));
575 MyVector3d d2Xdxideta(d2xyzdxideta[qp](0), d2xyzdxideta[qp](1), d2xyzdxideta[qp](2));
576
577
578 MyMatrix2d b;
579 b <<
580 n.dot(d2Xdxi2), n.dot(d2Xdxideta),
581 n.dot(d2Xdxideta), n.dot(d2Xdeta2);
582
583 MyVector3d dndxi = -b(0,0)*F0it.col(0) - b(0,1)*F0it.col(1);
584 MyVector3d dndeta = -b(1,0)*F0it.col(0) - b(1,1)*F0it.col(1);
585
586 MyMatrix2d bhat;
587 bhat <<
588 F0it.col(1).dot(dndeta), -F0it.col(0).dot(dndeta),
589 -F0it.col(1).dot(dndxi), F0it.col(0).dot(dndxi);
590
591 MyMatrix2d bc;
592 bc = bhat*C0;
593
594 // Mean curvature
595 Real H = 0.5*(dndxi.dot(F0it.col(0))+dndeta.dot(F0it.col(1)));
596
597 // Quadrature point reference coordinates
598 Real xi = qrule.qp(qp)(0);
599 Real eta = qrule.qp(qp)(1);
600
601 // Preassemble the MITC4 shear strain matrix for all nodes as they involve
602 // cross references to midside nodes.
603 // The QUAD4 element has nodes X1,X2,X3,X4 with coordinates (xi,eta)
604 // in the reference element: (-1,-1),(1,-1),(1,1),(-1,1).
605 // The midside nodes are denoted A1=(X1+X2)/2, B2=(X2+X3)/2, A2=(X3+X4)/2, B1=(X4+X1)/2.
606
607 // Normals at the midside nodes (average of normals at the edge corners).
608 // Multiplication by the assumed shear strain shape function.
609 MyVector3d nA1 = 0.5*(Qnode[0].col(2)+Qnode[1].col(2));
610 nA1 /= nA1.norm();
611 nA1 *= (1-eta)/4;
612 MyVector3d nB2 = 0.5*(Qnode[1].col(2)+Qnode[2].col(2));
613 nB2 /= nB2.norm();
614 nB2 *= (1+xi)/4;
615 MyVector3d nA2 = 0.5*(Qnode[2].col(2)+Qnode[3].col(2));
616 nA2 /= nA2.norm();
617 nA2 *= (1+eta)/4;
618 MyVector3d nB1 = 0.5*(Qnode[3].col(2)+Qnode[0].col(2));
619 nB1 /= nB1.norm();
620 nB1 *= (1-xi)/4;
621
622 // Edge tangents
623 MyVector3d aA1 = 0.5*(X2-X1);
624 MyVector3d aA2 = 0.5*(X3-X4);
625 MyVector3d aB1 = 0.5*(X4-X1);
626 MyVector3d aB2 = 0.5*(X3-X2);
627
628 // Contribution of the rotational dofs to the shear strain
629 MyVector2d AS1A1(-aA1.dot(Qnode[0].col(1)), aA1.dot(Qnode[0].col(0)));
630 MyVector2d AS2A1(-aA1.dot(Qnode[1].col(1)), aA1.dot(Qnode[1].col(0)));
631 AS1A1 *= (1-eta)/4;
632 AS2A1 *= (1-eta)/4;
633
634 MyVector2d AS1A2(-aA2.dot(Qnode[3].col(1)), aA2.dot(Qnode[3].col(0)));
635 MyVector2d AS2A2(-aA2.dot(Qnode[2].col(1)), aA2.dot(Qnode[2].col(0)));
636 AS1A2 *= (1+eta)/4;
637 AS2A2 *= (1+eta)/4;
638
639 MyVector2d AS1B1(-aB1.dot(Qnode[0].col(1)), aB1.dot(Qnode[0].col(0)));
640 MyVector2d AS2B1(-aB1.dot(Qnode[3].col(1)), aB1.dot(Qnode[3].col(0)));
641 AS1B1 *= (1-xi)/4;
642 AS2B1 *= (1-xi)/4;
643
644 MyVector2d AS1B2(-aB2.dot(Qnode[1].col(1)), aB2.dot(Qnode[1].col(0)));
645 MyVector2d AS2B2(-aB2.dot(Qnode[2].col(1)), aB2.dot(Qnode[2].col(0)));
646 AS1B2 *= (1+xi)/4;
647 AS2B2 *= (1+xi)/4;
648
649 // Store previous quantities in the shear strain matrices for each node
650 std::vector<MyMatrixXd> Bcnode;
651 MyMatrixXd Bc(2, 5);
652 // Node 1
653 Bc.block<1,3>(0,0) = -nA1.transpose();
654 Bc.block<1,2>(0,3) = AS1A1.transpose();
655 Bc.block<1,3>(1,0) = -nB1.transpose();
656 Bc.block<1,2>(1,3) = AS1B1.transpose();
657 Bcnode.push_back(Bc);
658 // Node 2
659 Bc.block<1,3>(0,0) = nA1.transpose();
660 Bc.block<1,2>(0,3) = AS2A1.transpose();
661 Bc.block<1,3>(1,0) = -nB2.transpose();
662 Bc.block<1,2>(1,3) = AS1B2.transpose();
663 Bcnode.push_back(Bc);
664 // Node 3
665 Bc.block<1,3>(0,0) = nA2.transpose();
666 Bc.block<1,2>(0,3) = AS2A2.transpose();
667 Bc.block<1,3>(1,0) = nB2.transpose();
668 Bc.block<1,2>(1,3) = AS2B2.transpose();
669 Bcnode.push_back(Bc);
670 // Node 4
671 Bc.block<1,3>(0,0) = -nA2.transpose();
672 Bc.block<1,2>(0,3) = AS1A2.transpose();
673 Bc.block<1,3>(1,0) = nB1.transpose();
674 Bc.block<1,2>(1,3) = AS2B1.transpose();
675 Bcnode.push_back(Bc);
676
677 // Loop over all pairs of nodes I,J.
678 for (unsigned int i=0; i<n_var_dofs; ++i)
679 {
680 // Matrix B0, zeroth order (through thickness) membrane-bending strain
681 Real C1i = dphidxi[i][qp]*C0(0,0) + dphideta[i][qp]*C0(1,0);
682 Real C2i = dphidxi[i][qp]*C0(0,1) + dphideta[i][qp]*C0(1,1);
683
684 MyMatrixXd B0I(3, 5);
685 B0I = MyMatrixXd::Zero(3, 5);
686 B0I.block<1,3>(0,0) = C1i*Q.col(0).transpose();
687 B0I.block<1,3>(1,0) = C2i*Q.col(1).transpose();
688 B0I.block<1,3>(2,0) = C2i*Q.col(0).transpose()+C1i*Q.col(1).transpose();
689
690 // Matrix B1, first order membrane-bending strain
691 Real bc1i = dphidxi[i][qp]*bc(0,0) + dphideta[i][qp]*bc(1,0);
692 Real bc2i = dphidxi[i][qp]*bc(0,1) + dphideta[i][qp]*bc(1,1);
693
694 MyVector2d V1i(-Q.col(0).dot(Qnode[i].col(1)),
695 Q.col(0).dot(Qnode[i].col(0)));
696
697 MyVector2d V2i(-Q.col(1).dot(Qnode[i].col(1)),
698 Q.col(1).dot(Qnode[i].col(0)));
699
700 MyMatrixXd B1I(3,5);
701 B1I = MyMatrixXd::Zero(3,5);
702 B1I.block<1,3>(0,0) = bc1i*Q.col(0).transpose();
703 B1I.block<1,3>(1,0) = bc2i*Q.col(1).transpose();
704 B1I.block<1,3>(2,0) = bc2i*Q.col(0).transpose()+bc1i*Q.col(1).transpose();
705
706 B1I.block<1,2>(0,3) = C1i*V1i.transpose();
707 B1I.block<1,2>(1,3) = C2i*V2i.transpose();
708 B1I.block<1,2>(2,3) = C2i*V1i.transpose()+C1i*V2i.transpose();
709
710 // Matrix B2, second order membrane-bending strain
711 MyMatrixXd B2I(3,5);
712 B2I = MyMatrixXd::Zero(3,5);
713
714 B2I.block<1,2>(0,3) = bc1i*V1i.transpose();
715 B2I.block<1,2>(1,3) = bc2i*V2i.transpose();
716 B2I.block<1,2>(2,3) = bc2i*V1i.transpose()+bc1i*V2i.transpose();
717
718 // Matrix Bc0, zeroth order shear strain
719 MyMatrixXd Bc0I(2,5);
720 Bc0I = C0.transpose()*Bcnode[i];
721
722 // Matrix Bc1, first order shear strain
723 MyMatrixXd Bc1I(2,5);
724 Bc1I = bc.transpose()*Bcnode[i];
725
726 // Drilling dof (in-plane rotation)
727 MyVector2d BdxiI(dphidxi[i][qp],dphideta[i][qp]);
728 MyVector2d BdI = C0.transpose()*BdxiI;
729
730 for (unsigned int j=0; j<n_var_dofs; ++j)
731 {
732
733 // Matrix B0, zeroth order membrane-bending strain
734 Real C1j = dphidxi[j][qp]*C0(0,0) + dphideta[j][qp]*C0(1,0);
735 Real C2j = dphidxi[j][qp]*C0(0,1) + dphideta[j][qp]*C0(1,1);
736
737 MyMatrixXd B0J(3,5);
738 B0J = MyMatrixXd::Zero(3,5);
739 B0J.block<1,3>(0,0) = C1j*Q.col(0).transpose();
740 B0J.block<1,3>(1,0) = C2j*Q.col(1).transpose();
741 B0J.block<1,3>(2,0) = C2j*Q.col(0).transpose()+C1j*Q.col(1).transpose();
742
743 // Matrix B1, first order membrane-bending strain
744 Real bc1j = dphidxi[j][qp]*bc(0,0) + dphideta[j][qp]*bc(1,0);
745 Real bc2j = dphidxi[j][qp]*bc(0,1) + dphideta[j][qp]*bc(1,1);
746
747 MyVector2d V1j(-Q.col(0).dot(Qnode[j].col(1)),
748 Q.col(0).dot(Qnode[j].col(0)));
749
750 MyVector2d V2j(-Q.col(1).dot(Qnode[j].col(1)),
751 Q.col(1).dot(Qnode[j].col(0)));
752
753 MyMatrixXd B1J(3,5);
754 B1J = MyMatrixXd::Zero(3,5);
755 B1J.block<1,3>(0,0) = bc1j*Q.col(0).transpose();
756 B1J.block<1,3>(1,0) = bc2j*Q.col(1).transpose();
757 B1J.block<1,3>(2,0) = bc2j*Q.col(0).transpose()+bc1j*Q.col(1).transpose();
758
759 B1J.block<1,2>(0,3) = C1j*V1j.transpose();
760 B1J.block<1,2>(1,3) = C2j*V2j.transpose();
761 B1J.block<1,2>(2,3) = C2j*V1j.transpose()+C1j*V2j.transpose();
762
763 // Matrix B2, second order membrane-bending strain
764 MyMatrixXd B2J(3,5);
765 B2J = MyMatrixXd::Zero(3,5);
766
767 B2J.block<1,2>(0,3) = bc1j*V1j.transpose();
768 B2J.block<1,2>(1,3) = bc2j*V2j.transpose();
769 B2J.block<1,2>(2,3) = bc2j*V1j.transpose()+bc1j*V2j.transpose();
770
771 // Matrix Bc0, zeroth order shear strain
772 MyMatrixXd Bc0J(2, 5);
773 Bc0J = C0.transpose()*Bcnode[j];
774
775 // Matrix Bc1, first order shear strain
776 MyMatrixXd Bc1J(2, 5);
777 Bc1J = bc.transpose()*Bcnode[j];
778
779 // Drilling dof
780 MyVector2d BdxiJ(dphidxi[j][qp], dphideta[j][qp]);
781 MyVector2d BdJ = C0.transpose()*BdxiJ;
782
783 // The total stiffness matrix coupling the nodes
784 // I and J is a sum of membrane, bending and shear contributions.
785 MyMatrixXd local_KIJ(5, 5);
786 local_KIJ = JxW[qp] * (
787 B0I.transpose() * Hm * B0J
788 + B2I.transpose() * Hf * B0J
789 + B0I.transpose() * Hf * B2J
790 + B1I.transpose() * Hf * B1J
791 + 2*H * B0I.transpose() * Hf * B1J
792 + 2*H * B1I.transpose() * Hf * B0J
793 + Bc0I.transpose() * Hc0 * Bc0J
794 + Bc1I.transpose() * Hc1 * Bc1J
795 + 2*H * Bc0I.transpose() * Hc1 * Bc1J
796 + 2*H * Bc1I.transpose() * Hc1 * Bc0J
797 );
798
799 // Going from 5 to 6 dofs to add drilling dof
800 MyMatrixXd full_local_KIJ(6, 6);
801 full_local_KIJ = MyMatrixXd::Zero(6, 6);
802 full_local_KIJ.block<5,5>(0,0)=local_KIJ;
803
804 // Drilling dof stiffness contribution
805 // Note that in the original book, there is a coefficient of
806 // alpha between 1e-4 and 1e-7 to make the fictitious
807 // drilling stiffness small while preventing the stiffness
808 // matrix from being singular. For this problem, we can use
809 // alpha = 1 to also get a good result.
810 //
811 // The explicit conversion to Real here is to work
812 // around an Eigen+boost::float128 incompatibility.
813 full_local_KIJ(5,5) = Real(Hf(0,0)*JxW[qp]*BdI.transpose()*BdJ);
814
815 // Transform the stiffness matrix to global coordinates
816 MyMatrixXd global_KIJ(6,6);
817 MyMatrixXd TI(6,6);
818 TI = MyMatrixXd::Identity(6,6);
819 TI.block<3,3>(3,3) = Qnode[i].transpose();
820 MyMatrixXd TJ(6,6);
821 TJ = MyMatrixXd::Identity(6,6);
822 TJ.block<3,3>(3,3) = Qnode[j].transpose();
823 global_KIJ = TI.transpose()*full_local_KIJ*TJ;
824
825 // Insert the components of the coupling stiffness
826 // matrix KIJ into the corresponding directional
827 // submatrices.
828 for (unsigned int k=0;k<6;k++)
829 for (unsigned int l=0;l<6;l++)
830 Ke_var[k][l](i,j) += global_KIJ(k,l);
831 }
832 }
833
834 } // end of the quadrature point qp-loop
835
836 if (distributed_load)
837 {
838 // Loop on shell faces
839 for (unsigned int shellface=0; shellface<2; shellface++)
840 {
841 std::vector<boundary_id_type> bids;
842 mesh.get_boundary_info().shellface_boundary_ids(elem, shellface, bids);
843
844 for (std::size_t k=0; k<bids.size(); k++)
845 if (bids[k]==11) // sideset id for surface load
846 for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
847 for (unsigned int i=0; i<n_var_dofs; ++i)
848 Fe_w(i) -= JxW[qp] * phi[i][qp];
849 }
850 }
851
852 // The element matrix is now built for this element.
853 // Add it to the global matrix.
854
855 dof_map.constrain_element_matrix_and_vector (Ke,Fe,dof_indices);
856
857 system.matrix->add_matrix (Ke, dof_indices);
858 system.rhs->add_vector (Fe, dof_indices);
859
860 }
861
862 if (!distributed_load)
863 {
864 //Adding point load to the RHS
865
866 //Pinch position
867 Point C(0, 3, 3);
868
869 //Finish assembling rhs so we can set one value
870 system.rhs->close();
871
872 for (const auto & node : mesh.node_ptr_range())
873 if (((*node) - C).norm() < 1e-3)
874 system.rhs->set(node->dof_number(0, 2, 0), -q/4);
875 }
876
877 #else
878 // Avoid compiler warnings
879 libmesh_ignore(es, system_name);
880 #endif // defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
881 }
882