libMesh
miscellaneous_ex13.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 13 - Quad8 Shell Elements</h1>
19 // \author Sylvain Vallaghe
20 // \date 2017
21 //
22 // This example implements shell elements using 8-noded serendipity quadrilaterals and reduced integration.
23 // Full integration Quad8 shell elements are prone to locking issues.
24 // The reduced integration version "eliminates locking in most situations
25 // although it introduces two spureous mechanisms. Fortunately,
26 // these mechanisms disappear after assembly of the global stiffness matrix and the element
27 // can be considered safe for practical purposes".
28 // (source: Onate, Structural Analysis with the Finite Element Method).
29 // The "pinched cylinder" problem is solved and the solution
30 // is compared to analytical values at selected points.
31 
32 // C++ include files that we need
33 #include <iostream>
34 
35 // LibMesh includes
36 #include "libmesh/libmesh.h"
37 #include "libmesh/mesh.h"
38 #include "libmesh/linear_implicit_system.h"
39 #include "libmesh/equation_systems.h"
40 #include "libmesh/fe.h"
41 #include "libmesh/quadrature.h"
42 #include "libmesh/node.h"
43 #include "libmesh/elem.h"
44 #include "libmesh/dof_map.h"
45 #include "libmesh/quadrature_gauss.h"
46 #include "libmesh/vector_value.h"
47 #include "libmesh/tensor_value.h"
48 #include "libmesh/dense_matrix.h"
49 #include "libmesh/dense_submatrix.h"
50 #include "libmesh/dense_vector.h"
51 #include "libmesh/dense_subvector.h"
52 #include "libmesh/sparse_matrix.h"
53 #include "libmesh/numeric_vector.h"
54 #include "libmesh/exodusII_io.h"
55 #include "libmesh/dirichlet_boundaries.h"
56 #include "libmesh/zero_function.h"
57 #include "libmesh/linear_solver.h"
58 #include "libmesh/getpot.h"
59 #include "libmesh/enum_solver_package.h"
60 #include "libmesh/enum_solver_type.h"
61 #include "libmesh/parallel.h"
62 
63 // Eigen includes
64 #ifdef LIBMESH_HAVE_EIGEN
65 #include "libmesh/ignore_warnings.h"
66 # include <Eigen/Dense>
67 #include "libmesh/restore_warnings.h"
68 
69 // Use Eigen dense numerics with libMesh::Real
70 typedef Eigen::Matrix<libMesh::Real, 2, 2> MyMatrix2d;
71 typedef Eigen::Matrix<libMesh::Real, 3, 3> MyMatrix3d;
72 typedef Eigen::Matrix<libMesh::Real, 2, 1> MyVector2d;
73 typedef Eigen::Matrix<libMesh::Real, 3, 1> MyVector3d;
74 typedef Eigen::Matrix<libMesh::Real, Eigen::Dynamic, Eigen::Dynamic> MyMatrixXd;
75 #endif
76 
77 // Bring in everything from the libMesh namespace
78 using namespace libMesh;
79 
80 // Function prototype. This is the function that will assemble
81 // the stiffness matrix and the right-hand-side vector ready
82 // for solution.
84  const std::string & system_name);
85 
86 // Begin the main program.
87 int main (int argc, char ** argv)
88 {
89  // Initialize libMesh.
90  LibMeshInit init (argc, argv);
91 
92  // Skip this 3D example if libMesh was compiled as 1D/2D-only.
93  libmesh_example_requires (3 == LIBMESH_DIM, "3D support");
94 
95  // We use Dirichlet boundary conditions here
96 #ifndef LIBMESH_ENABLE_DIRICHLET
97  libmesh_example_requires(false, "--enable-dirichlet");
98 #endif
99 
100  // Our input mesh here is in ExodusII format
101 #ifndef LIBMESH_HAVE_EXODUS_API
102  libmesh_example_requires (false, "ExodusII support");
103 #endif
104 
105 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
106  libmesh_example_requires (false, "second derivatives enabled");
107 #endif
108 
109  // This example does a bunch of linear algebra during assembly, and
110  // therefore requires Eigen.
111 #ifndef LIBMESH_HAVE_EIGEN
112  libmesh_example_requires(false, "--enable-eigen");
113 #endif
114 
115  // This example converts between ExodusII and XDR files, therefore
116  // it requires XDR support in libmesh.
117 #ifndef LIBMESH_HAVE_XDR
118  libmesh_example_requires (false, "XDR support");
119 #endif
120 
121  // This examples requires parallel minloc() support, which we don't
122  // implement yet for float128
123 #ifdef LIBMESH_DEFAULT_QUADRUPLE_PRECISION
124  libmesh_example_requires (false, "--disable-quadruple-precision");
125 #endif
126 
127  // Read the "distributed_load" flag from the command line
128  const int distributed_load = libMesh::command_line_next("-distributed_load", 0);
129 
130  {
131  Mesh mesh (init.comm(), 3);
132 
133  // To confirm that both ExodusII and Xdr formats work for shell
134  // meshes, we read in cylinder.exo, then write out cylinder.xdr,
135  // then read in cylinder.exo again below and use that for the rest
136  // of the example.
137  mesh.read("cylinder.exo");
138  mesh.write("cylinder.xdr");
139  }
140  Mesh mesh (init.comm(), 3);
141  mesh.read("cylinder.xdr");
142 
143  // Print information about the mesh to the screen.
144  mesh.print_info();
145 
146  // Create an equation systems object.
147  EquationSystems equation_systems (mesh);
148 
149  // Declare the system and its variables.
150  // Create a linear implicit system named "Shell".
151  LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem> ("Shell");
152 
153  // Add the three displacement variables "u", "v", "w",
154  // and the three rotational variables "theta_x", "theta_y", "theta_z".
155  // All variables are second order.
156  system.add_variable ("u",SECOND,LAGRANGE);
157  system.add_variable ("v",SECOND,LAGRANGE);
158  system.add_variable ("w",SECOND,LAGRANGE);
159  system.add_variable ("theta_x",SECOND,LAGRANGE);
160  system.add_variable ("theta_y",SECOND,LAGRANGE);
161  system.add_variable ("theta_z",SECOND,LAGRANGE);
162 
163  // Give the system a pointer to the matrix and rhs assembly
164  // function.
166 
167  // Use the parameters of the equation systems object to
168  // tell the shell system about the material properties, the
169  // shell thickness, and the external load.
170  const Real h = 0.03;
171  const Real E = 3e10;
172  const Real nu = 0.3;
173  const Real q = 1;
174  equation_systems.parameters.set<Real> ("thickness") = h;
175  equation_systems.parameters.set<Real> ("young's modulus") = E;
176  equation_systems.parameters.set<Real> ("poisson ratio") = nu;
177  equation_systems.parameters.set<Real> ("point load") = q;
178  equation_systems.parameters.set<bool>("distributed load") = (distributed_load != 0);
179 
180  // Dirichlet conditions for the pinched cylinder problem.
181  // Only one 8th of the cylinder is considered using symmetry considerations.
182  // The cylinder longitudinal axis is the y-axis.
183  // The four corners of the surface are named A(3,0,0), B(3,3,0), C(0,3,3), D(0,0,3).
184  // The point load (pinch) is applied at C in the -z direction.
185  // Edge AD is the actual edge of the cylinder and is rigid in the xz-plane.
186  // Other edges have symmetric boundary conditions.
187 
188 #ifdef LIBMESH_ENABLE_DIRICHLET
189  // AB w, theta_x, theta_y
190  {
191  ZeroFunction<> zf;
192 
193  // Most DirichletBoundary users will want to supply a "locally
194  // indexed" functor
195  DirichletBoundary dirichlet_bc
196  (/*boundary_ids =*/{7},/*variables =*/{2,3,4}, zf,
198  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
199  }
200  // BC v, theta_x, theta_z
201  {
202  ZeroFunction<> zf;
203 
204  DirichletBoundary dirichlet_bc
205  (/*boundary_ids =*/{8}, /*variables =*/{1,3,5}, zf,
207  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
208  }
209  // CD u, theta_y, theta_z
210  {
211  ZeroFunction<> zf;
212 
213  DirichletBoundary dirichlet_bc
214  (/*boundary_ids =*/{9}, /*variables =*/{0,4,5}, zf,
216  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
217  }
218  // AD u, w, theta_y
219  {
220  ZeroFunction<> zf;
221 
222  DirichletBoundary dirichlet_bc
223  (/*boundary_ids =*/{10}, /*variables =*/{0,2,4}, zf,
225  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
226  }
227 #endif // LIBMESH_ENABLE_DIRICHLET
228 
229  // Initialize the data structures for the equation system.
230  equation_systems.init();
231 
232  // Print information about the system to the screen.
233  equation_systems.print_info();
234 
235  // Solve the linear system.
236  system.solve();
237 
238  // After solving the system, write the solution to an
239  // ExodusII output file ready for import in, e.g.,
240  // Paraview.
241  ExodusII_IO(mesh).write_equation_systems ("out.e", equation_systems);
242 
243  // Compare with analytical solution for point load
244  if (distributed_load==0)
245  {
246  // Find the node nearest point C.
247  Node * node_C = nullptr;
248  Point point_C(0, 3, 3);
249  {
250  Real nearest_dist_sq = std::numeric_limits<Real>::max();
251 
252  // Find the closest local node. On a DistributedMesh we may
253  // not even know about the existence of closer non-local
254  // nodes.
255  for (const auto & node : mesh.local_node_ptr_range())
256  {
257  const Real dist_sq = (*node - point_C).norm_sq();
258  if (dist_sq < nearest_dist_sq)
259  {
260  nearest_dist_sq = dist_sq;
261  node_C = node;
262  }
263  }
264 
265  // Check with other processors to see if any found a closer node
266  unsigned int minrank = 0;
267  system.comm().minloc(nearest_dist_sq, minrank);
268 
269  // Broadcast the ID of the closest node, so every processor can
270  // see for certain whether they have it or not.
271  dof_id_type nearest_node_id = 0;
272  if (system.processor_id() == minrank)
273  nearest_node_id = node_C->id();
274  system.comm().broadcast(nearest_node_id, minrank);
275  node_C = mesh.query_node_ptr(nearest_node_id);
276  }
277 
278  // Evaluate the z-displacement "w" at the node nearest C.
279  Number w = 0;
280 
281  // If we know about the closest node, and if we also own the DoFs
282  // on that node, then we can evaluate the solution at that node.
283  if (node_C)
284  {
285  const unsigned int w_var = system.variable_number ("w");
286  dof_id_type w_dof = node_C->dof_number (system.number(), w_var, 0);
287  if (w_dof >= system.get_dof_map().first_dof() &&
288  w_dof < system.get_dof_map().end_dof())
289  w = system.current_solution(w_dof);
290  }
291  system.comm().sum(w);
292 
293 
294  Number w_C_bar = -E*h*w/q;
295  const Real w_C_bar_analytic = 164.24;
296 
297  // Print the finite element solution and the analytic
298  // prediction to the screen.
299  libMesh::out << "z-displacement of the point C: " << w_C_bar << std::endl;
300  libMesh::out << "Analytic solution: " << w_C_bar_analytic << std::endl;
301 
302  // Evaluate the y-displacement "v" at point D. This time we'll
303  // evaluate at the exact point, not just the closest node.
304  Point point_D(0, 0, 3);
305  const unsigned int v_var = system.variable_number ("v");
306  Number v = system.point_value(v_var, point_D);
307 
308  Number v_D_bar = E*h*v/q;
309  const Real v_D_bar_analytic = 4.114;
310 
311  // Print the finite element solution and the analytic
312  // prediction to the screen.
313  libMesh::out << "y-displacement of the point D: " << v_D_bar << std::endl;
314  libMesh::out << "Analytic solution: " << v_D_bar_analytic << std::endl;
315  }
316 
317  // All done.
318  return 0;
319 }
320 
321 
322 
323 // We now define the matrix and rhs vector assembly function
324 // for the shell system.
326  const std::string & system_name)
327 {
328  // This example requires Eigen to actually work, but we should still
329  // let it compile and throw a runtime error if you don't.
330 #if defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
331  // It is a good idea to make sure we are assembling
332  // the proper system.
333  libmesh_assert_equal_to (system_name, "Shell");
334 
335  // Get a constant reference to the mesh object.
336  const MeshBase & mesh = es.get_mesh();
337  const unsigned int dim = mesh.mesh_dimension();
338 
339  // Get a reference to the shell system object.
340  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> (system_name);
341 
342  // Get the shell parameters that we need during assembly.
343  const Real h = es.parameters.get<Real> ("thickness");
344  const Real E = es.parameters.get<Real> ("young's modulus");
345  const Real nu = es.parameters.get<Real> ("poisson ratio");
346  const Real q = es.parameters.get<Real> ("point load");
347  const bool distributed_load = es.parameters.get<bool> ("distributed load");
348 
349  // The membrane elastic matrix.
350  MyMatrix3d Hm;
351  Hm <<
352  1., nu, 0.,
353  nu, 1., 0.,
354  0., 0., 0.5 * (1-nu);
355  Hm *= h * E/(1-nu*nu);
356 
357  // The bending elastic matrix.
358  MyMatrix3d Hf;
359  Hf <<
360  1., nu, 0.,
361  nu, 1., 0.,
362  0., 0., 0.5 * (1-nu);
363  Hf *= h*h*h/12 * E/(1-nu*nu);
364 
365  // The shear elastic matrices.
366  MyMatrix2d Hc0 = MyMatrix2d::Identity();
367  Hc0 *= h * 5./6*E/(2*(1+nu));
368 
369  MyMatrix2d Hc1 = MyMatrix2d::Identity();
370  Hc1 *= h*h*h/12 * 5./6*E/(2*(1+nu));
371 
372  // Get the Finite Element type, this will be
373  // the same for all variables.
374  FEType fe_type = system.variable_type (0);
375 
376  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
377  QGauss qrule (dim, SECOND); // Reduced integration, 2x2 gauss points (instead of 3x3 for full integration)
378  fe->attach_quadrature_rule (&qrule);
379 
380  // The element Jacobian * quadrature weight at each integration point.
381  const std::vector<Real> & JxW = fe->get_JxW();
382 
383  // The element shape function and its derivatives evaluated at the
384  // quadrature points.
385  const std::vector<RealGradient> & dxyzdxi = fe->get_dxyzdxi();
386  const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();
387  const std::vector<RealGradient> & d2xyzdxi2 = fe->get_d2xyzdxi2();
388  const std::vector<RealGradient> & d2xyzdeta2 = fe->get_d2xyzdeta2();
389  const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();
390  const std::vector<std::vector<Real>> & dphidxi = fe->get_dphidxi();
391  const std::vector<std::vector<Real>> & dphideta = fe->get_dphideta();
392  const std::vector<std::vector<Real>> & phi = fe->get_phi();
393 
394  // A reference to the DofMap object for this system. The DofMap
395  // object handles the index translation from node and element numbers
396  // to degree of freedom numbers.
397  const DofMap & dof_map = system.get_dof_map();
398 
399  // The global system matrix
400  SparseMatrix<Number> & matrix = system.get_system_matrix();
401 
402  // Define data structures to contain the element stiffness matrix.
404  DenseSubMatrix<Number> Ke_var[6][6] =
405  {
418  };
419 
420  // Define data structures to contain the element rhs vector.
422  DenseSubVector<Number> Fe_w(Fe);
423 
424  std::vector<dof_id_type> dof_indices;
425  std::vector<std::vector<dof_id_type>> dof_indices_var(6);
426 
427  // Now we will loop over all the elements in the mesh. We will
428  // compute the element matrix and right-hand-side contribution.
429  for (const auto & elem : mesh.active_local_element_ptr_range())
430  {
431  dof_map.dof_indices (elem, dof_indices);
432  for (unsigned int var=0; var<6; var++)
433  dof_map.dof_indices (elem, dof_indices_var[var], var);
434 
435  const unsigned int n_dofs = dof_indices.size();
436  const unsigned int n_var_dofs = dof_indices_var[0].size();
437 
438  // First compute element data at the nodes
439  std::vector<Point> nodes;
440  for (unsigned int i=0; i<elem->n_nodes(); ++i)
441  nodes.push_back(elem->reference_elem()->node_ref(i));
442  fe->reinit (elem, &nodes);
443 
444  //Store local orthonormal basis at the nodes
445  std::vector<MyMatrix3d> Qnode;
446  for (unsigned int i=0; i<elem->n_nodes(); ++i)
447  {
448  MyVector3d a1;
449  a1 << dxyzdxi[i](0), dxyzdxi[i](1), dxyzdxi[i](2);
450  MyVector3d a2;
451  a2 << dxyzdeta[i](0), dxyzdeta[i](1), dxyzdeta[i](2);
452  MyVector3d n;
453  n = a1.cross(a2);
454  n /= n.norm();
455 
456  Real nx = n(0);
457  Real ny = n(1);
458  Real C = n(2);
459  if (std::abs(1.+C)<1e-6)
460  {
461  MyMatrix3d Q;
462  Q <<
463  1, 0, 0,
464  0, -1, 0,
465  0, 0, -1;
466  Qnode.push_back(Q);
467  }
468  else
469  {
470  MyMatrix3d Q;
471  Q <<
472  C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
473  -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
474  -nx, -ny, C;
475  Qnode.push_back(Q);
476  }
477  }
478 
479  Ke.resize (n_dofs, n_dofs);
480  for (unsigned int var_i=0; var_i<6; var_i++)
481  for (unsigned int var_j=0; var_j<6; var_j++)
482  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
483 
484  Fe.resize(n_dofs);
485  Fe_w.reposition(2*n_var_dofs,n_var_dofs);
486 
487  // Reinit element data at the regular Gauss quadrature points
488  fe->reinit (elem);
489 
490  // Now we will build the element matrix and right-hand-side.
491  for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
492  {
493 
494  //Covariant basis at the quadrature point
495  MyVector3d a1;
496  a1 << dxyzdxi[qp](0), dxyzdxi[qp](1), dxyzdxi[qp](2);
497  MyVector3d a2;
498  a2 << dxyzdeta[qp](0), dxyzdeta[qp](1), dxyzdeta[qp](2);
499  MyVector3d n;
500  n = a1.cross(a2);
501  n /= n.norm();
502  MyMatrix3d F0;
503  F0 <<
504  a1(0), a2(0), n(0),
505  a1(1), a2(1), n(1),
506  a1(2), a2(2), n(2);
507 
508  //Contravariant basis
509  MyMatrix3d F0it;
510  F0it = F0.inverse().transpose();
511 
512  //Local orthonormal basis at the quadrature point
513  Real nx = n(0);
514  Real ny = n(1);
515  Real C = n(2);
516  MyMatrix3d Q;
517  if (std::abs(1.+C) < 1e-6)
518  {
519  Q <<
520  1, 0, 0,
521  0, -1, 0,
522  0, 0, -1;
523  }
524  else
525  {
526  Q <<
527  C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
528  -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
529  -nx, -ny, C;
530  }
531 
532  MyMatrix2d C0;
533  C0 = F0it.block<3,2>(0,0).transpose()*Q.block<3,2>(0,0);
534 
535  // Normal derivatives in reference coordinates
536  MyVector3d d2Xdxi2(d2xyzdxi2[qp](0), d2xyzdxi2[qp](1), d2xyzdxi2[qp](2));
537  MyVector3d d2Xdeta2(d2xyzdeta2[qp](0), d2xyzdeta2[qp](1), d2xyzdeta2[qp](2));
538  MyVector3d d2Xdxideta(d2xyzdxideta[qp](0), d2xyzdxideta[qp](1), d2xyzdxideta[qp](2));
539 
540  MyMatrix2d b;
541  b <<
542  n.dot(d2Xdxi2), n.dot(d2Xdxideta),
543  n.dot(d2Xdxideta), n.dot(d2Xdeta2);
544 
545  MyVector3d dndxi = -b(0,0)*F0it.col(0) - b(0,1)*F0it.col(1);
546  MyVector3d dndeta = -b(1,0)*F0it.col(0) - b(1,1)*F0it.col(1);
547 
548  MyMatrix2d bhat;
549  bhat <<
550  F0it.col(1).dot(dndeta), -F0it.col(0).dot(dndeta),
551  -F0it.col(1).dot(dndxi), F0it.col(0).dot(dndxi);
552 
553  MyMatrix2d bc;
554  bc = bhat*C0;
555 
556  // Mean curvature
557  Real H = 0.5*(dndxi.dot(F0it.col(0))+dndeta.dot(F0it.col(1)));
558 
559  // Loop over all pairs of nodes I,J.
560  for (unsigned int i=0; i<n_var_dofs; ++i)
561  {
562  // Matrix B0, zeroth order (through thickness) membrane-bending strain
563  Real C1i = dphidxi[i][qp]*C0(0,0) + dphideta[i][qp]*C0(1,0);
564  Real C2i = dphidxi[i][qp]*C0(0,1) + dphideta[i][qp]*C0(1,1);
565 
566  MyMatrixXd B0I(3, 5);
567  B0I = MyMatrixXd::Zero(3, 5);
568  B0I.block<1,3>(0,0) = C1i*Q.col(0).transpose();
569  B0I.block<1,3>(1,0) = C2i*Q.col(1).transpose();
570  B0I.block<1,3>(2,0) = C2i*Q.col(0).transpose()+C1i*Q.col(1).transpose();
571 
572  // Matrix B1, first order membrane-bending strain
573  Real bc1i = dphidxi[i][qp]*bc(0,0) + dphideta[i][qp]*bc(1,0);
574  Real bc2i = dphidxi[i][qp]*bc(0,1) + dphideta[i][qp]*bc(1,1);
575 
576  MyVector2d V1i(-Q.col(0).dot(Qnode[i].col(1)),
577  Q.col(0).dot(Qnode[i].col(0)));
578 
579  MyVector2d V2i(-Q.col(1).dot(Qnode[i].col(1)),
580  Q.col(1).dot(Qnode[i].col(0)));
581 
582  MyMatrixXd B1I(3,5);
583  B1I = MyMatrixXd::Zero(3,5);
584  B1I.block<1,3>(0,0) = bc1i*Q.col(0).transpose();
585  B1I.block<1,3>(1,0) = bc2i*Q.col(1).transpose();
586  B1I.block<1,3>(2,0) = bc2i*Q.col(0).transpose()+bc1i*Q.col(1).transpose();
587 
588  B1I.block<1,2>(0,3) = C1i*V1i.transpose();
589  B1I.block<1,2>(1,3) = C2i*V2i.transpose();
590  B1I.block<1,2>(2,3) = C2i*V1i.transpose()+C1i*V2i.transpose();
591 
592  // Matrix B2, second order membrane-bending strain
593  MyMatrixXd B2I(3,5);
594  B2I = MyMatrixXd::Zero(3,5);
595 
596  B2I.block<1,2>(0,3) = bc1i*V1i.transpose();
597  B2I.block<1,2>(1,3) = bc2i*V2i.transpose();
598  B2I.block<1,2>(2,3) = bc2i*V1i.transpose()+bc1i*V2i.transpose();
599 
600  // Matrix Bc0, zeroth order shear strain
601  MyMatrixXd Bc0I(2,5);
602  Bc0I = MyMatrixXd::Zero(2,5);
603  Bc0I.block<1,3>(0,0) = C1i*Q.col(2).transpose();
604  Bc0I.block<1,3>(1,0) = C2i*Q.col(2).transpose();
605  Bc0I.block<1,2>(0,3) = phi[i][qp]*V1i.transpose();
606  Bc0I.block<1,2>(1,3) = phi[i][qp]*V2i.transpose();
607 
608  // Matrix Bc1, first order shear strain
609  MyMatrixXd Bc1I(2,5);
610  Bc1I = MyMatrixXd::Zero(2,5);
611  Bc1I.block<1,3>(0,0) = bc1i*Q.col(2).transpose();
612  Bc1I.block<1,3>(1,0) = bc2i*Q.col(2).transpose();
613 
614  // Drilling dof (in-plane rotation)
615  MyVector2d BdxiI(dphidxi[i][qp],dphideta[i][qp]);
616  MyVector2d BdI = C0.transpose()*BdxiI;
617 
618  for (unsigned int j=0; j<n_var_dofs; ++j)
619  {
620 
621  // Matrix B0, zeroth order membrane-bending strain
622  Real C1j = dphidxi[j][qp]*C0(0,0) + dphideta[j][qp]*C0(1,0);
623  Real C2j = dphidxi[j][qp]*C0(0,1) + dphideta[j][qp]*C0(1,1);
624 
625  MyMatrixXd B0J(3,5);
626  B0J = MyMatrixXd::Zero(3,5);
627  B0J.block<1,3>(0,0) = C1j*Q.col(0).transpose();
628  B0J.block<1,3>(1,0) = C2j*Q.col(1).transpose();
629  B0J.block<1,3>(2,0) = C2j*Q.col(0).transpose()+C1j*Q.col(1).transpose();
630 
631  // Matrix B1, first order membrane-bending strain
632  Real bc1j = dphidxi[j][qp]*bc(0,0) + dphideta[j][qp]*bc(1,0);
633  Real bc2j = dphidxi[j][qp]*bc(0,1) + dphideta[j][qp]*bc(1,1);
634 
635  MyVector2d V1j(-Q.col(0).dot(Qnode[j].col(1)),
636  Q.col(0).dot(Qnode[j].col(0)));
637 
638  MyVector2d V2j(-Q.col(1).dot(Qnode[j].col(1)),
639  Q.col(1).dot(Qnode[j].col(0)));
640 
641  MyMatrixXd B1J(3,5);
642  B1J = MyMatrixXd::Zero(3,5);
643  B1J.block<1,3>(0,0) = bc1j*Q.col(0).transpose();
644  B1J.block<1,3>(1,0) = bc2j*Q.col(1).transpose();
645  B1J.block<1,3>(2,0) = bc2j*Q.col(0).transpose()+bc1j*Q.col(1).transpose();
646 
647  B1J.block<1,2>(0,3) = C1j*V1j.transpose();
648  B1J.block<1,2>(1,3) = C2j*V2j.transpose();
649  B1J.block<1,2>(2,3) = C2j*V1j.transpose()+C1j*V2j.transpose();
650 
651  // Matrix B2, second order membrane-bending strain
652  MyMatrixXd B2J(3,5);
653  B2J = MyMatrixXd::Zero(3,5);
654 
655  B2J.block<1,2>(0,3) = bc1j*V1j.transpose();
656  B2J.block<1,2>(1,3) = bc2j*V2j.transpose();
657  B2J.block<1,2>(2,3) = bc2j*V1j.transpose()+bc1j*V2j.transpose();
658 
659  // Matrix Bc0, zeroth order shear strain
660  MyMatrixXd Bc0J(2, 5);
661  Bc0J = MyMatrixXd::Zero(2,5);
662  Bc0J.block<1,3>(0,0) = C1j*Q.col(2).transpose();
663  Bc0J.block<1,3>(1,0) = C2j*Q.col(2).transpose();
664  Bc0J.block<1,2>(0,3) = phi[j][qp]*V1j.transpose();
665  Bc0J.block<1,2>(1,3) = phi[j][qp]*V2j.transpose();
666 
667  // Matrix Bc1, first order shear strain
668  MyMatrixXd Bc1J(2, 5);
669  Bc1J = MyMatrixXd::Zero(2,5);
670  Bc1J.block<1,3>(0,0) = bc1j*Q.col(2).transpose();
671  Bc1J.block<1,3>(1,0) = bc2j*Q.col(2).transpose();
672 
673  // Drilling dof
674  MyVector2d BdxiJ(dphidxi[j][qp], dphideta[j][qp]);
675  MyVector2d BdJ = C0.transpose()*BdxiJ;
676 
677  // The total stiffness matrix coupling the nodes
678  // I and J is a sum of membrane, bending and shear contributions.
679  MyMatrixXd local_KIJ(5, 5);
680  local_KIJ = JxW[qp] * (
681  B0I.transpose() * Hm * B0J
682  + B2I.transpose() * Hf * B0J
683  + B0I.transpose() * Hf * B2J
684  + B1I.transpose() * Hf * B1J
685  + 2*H * B0I.transpose() * Hf * B1J
686  + 2*H * B1I.transpose() * Hf * B0J
687  + Bc0I.transpose() * Hc0 * Bc0J
688  + Bc1I.transpose() * Hc1 * Bc1J
689  + 2*H * Bc0I.transpose() * Hc1 * Bc1J
690  + 2*H * Bc1I.transpose() * Hc1 * Bc0J
691  );
692 
693  // Going from 5 to 6 dofs to add drilling dof
694  MyMatrixXd full_local_KIJ(6, 6);
695  full_local_KIJ = MyMatrixXd::Zero(6, 6);
696  full_local_KIJ.block<5,5>(0,0)=local_KIJ;
697 
698  // Drilling dof stiffness contribution
699  //
700  // The explicit conversion to Real here is to work
701  // around an Eigen+boost::float128 incompatibility.
702  full_local_KIJ(5,5) = Real(Hf(0,0)*JxW[qp]*BdI.transpose()*BdJ);
703 
704  // Transform the stiffness matrix to global coordinates
705  MyMatrixXd global_KIJ(6,6);
706  MyMatrixXd TI(6,6);
707  TI = MyMatrixXd::Identity(6,6);
708  TI.block<3,3>(3,3) = Qnode[i].transpose();
709  MyMatrixXd TJ(6,6);
710  TJ = MyMatrixXd::Identity(6,6);
711  TJ.block<3,3>(3,3) = Qnode[j].transpose();
712  global_KIJ = TI.transpose()*full_local_KIJ*TJ;
713 
714  // Insert the components of the coupling stiffness
715  // matrix KIJ into the corresponding directional
716  // submatrices.
717  for (unsigned int k=0;k<6;k++)
718  for (unsigned int l=0;l<6;l++)
719  Ke_var[k][l](i,j) += global_KIJ(k,l);
720  }
721  }
722 
723  } // end of the quadrature point qp-loop
724 
725  if (distributed_load)
726  {
727  // Loop on shell faces
728  for (unsigned int shellface=0; shellface<2; shellface++)
729  {
730  std::vector<boundary_id_type> bids;
731  mesh.get_boundary_info().shellface_boundary_ids(elem, shellface, bids);
732 
733  for (std::size_t k=0; k<bids.size(); k++)
734  if (bids[k]==11) // sideset id for surface load
735  for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
736  for (unsigned int i=0; i<n_var_dofs; ++i)
737  Fe_w(i) -= JxW[qp] * phi[i][qp];
738  }
739  }
740 
741  // The element matrix is now built for this element.
742  // Add it to the global matrix.
743 
744  dof_map.constrain_element_matrix_and_vector (Ke,Fe,dof_indices);
745 
746  matrix.add_matrix (Ke, dof_indices);
747  system.rhs->add_vector (Fe, dof_indices);
748 
749  }
750 
751  if (!distributed_load)
752  {
753  //Adding point load to the RHS
754 
755  //Pinch position
756  Point C(0, 3, 3);
757 
758  //Finish assembling rhs so we can set one value
759  system.rhs->close();
760 
761  for (const auto & node : mesh.node_ptr_range())
762  if (((*node) - C).norm() < 1e-3)
763  system.rhs->set(node->dof_number(0, 2, 0), -q/4);
764  }
765 
766 #else
767  // Avoid compiler warnings
768  libmesh_ignore(es, system_name);
769 #endif // defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
770 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
int main(int argc, char **argv)
T command_line_next(std::string name, T default_value)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:1011
This is the EquationSystems class.
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1025
A Node is like a Point, but with more information.
Definition: node.h:52
ConstFunction that simply returns 0.
Definition: zero_function.h:38
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:2254
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1992
void minloc(T &r, unsigned int &min_id) const
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
Eigen::Matrix< libMesh::Real, 3, 3 > MyMatrix3d
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
Eigen::Matrix< libMesh::Real, Eigen::Dynamic, Eigen::Dynamic > MyMatrixXd
void sum(T &r) const
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
MeshBase & mesh
Number point_value(unsigned int var, const Point &p, const bool insist_on_success=true, const NumericVector< Number > *sol=nullptr) const
Definition: system.C:2369
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
const Parallel::Communicator & comm() const
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:159
virtual void solve() override
Assembles & solves the linear system A*x=b.
const T_sys & get_system(std::string_view name) const
void assemble_shell(EquationSystems &es, const std::string &system_name)
Eigen::Matrix< libMesh::Real, 3, 3 > MyMatrix3d
This is the MeshBase class.
Definition: mesh_base.h:74
Eigen::Matrix< libMesh::Real, 2, 1 > MyVector2d
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:157
unsigned int variable_number(std::string_view var) const
Definition: system.C:1557
Defines a dense subvector for use in finite element computations.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:169
void reposition(const unsigned int ioff, const unsigned int n)
Changes the location of the subvector in the parent vector.
void libmesh_ignore(const Args &...)
unsigned int number() const
Definition: system.h:2269
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const T & get(std::string_view) const
Definition: parameters.h:426
dof_id_type id() const
Definition: dof_object.h:823
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Definition: mesh_base.C:1489
Defines a dense submatrix for use in Finite Element-type computations.
virtual const Node * query_node_ptr(const dof_id_type i) const =0
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1305
unsigned int n_points() const
Definition: quadrature.h:123
virtual void write(const std::string &name)=0
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:2109
Eigen::Matrix< libMesh::Real, 2, 2 > MyMatrix2d
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
Eigen::Matrix< libMesh::Real, Eigen::Dynamic, Eigen::Dynamic > MyMatrixXd
Definition: assembly.h:127
Eigen::Matrix< libMesh::Real, 3, 1 > MyVector3d
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2427
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
auto norm_sq(const T &a) -> decltype(std::norm(a))
Definition: tensor_tools.h:104
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
Eigen::Matrix< libMesh::Real, 2, 2 > MyMatrix2d
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:324
Parameters parameters
Data structure holding arbitrary parameters.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
virtual void init()
Initialize all the systems.
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:684
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:708
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
processor_id_type processor_id() const
Eigen::Matrix< libMesh::Real, 2, 1 > MyVector2d
const DofMap & get_dof_map() const
Definition: system.h:2293
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const SparseMatrix< Number > & get_system_matrix() const
Eigen::Matrix< libMesh::Real, 3, 1 > MyVector3d
uint8_t dof_id_type
Definition: id_types.h:67