libMesh
miscellaneous_ex13.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 - 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/libmesh_nullptr.h"
59 #include "libmesh/getpot.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 #endif
67 
68 // Bring in everything from the libMesh namespace
69 using namespace libMesh;
70 
71 // Function prototype. This is the function that will assemble
72 // the stiffness matrix and the right-hand-side vector ready
73 // for solution.
75  const std::string & system_name);
76 
77 // Begin the main program.
78 int main (int argc, char ** argv)
79 {
80  // Initialize libMesh.
81  LibMeshInit init (argc, argv);
82 
83  // Skip this 3D example if libMesh was compiled as 1D/2D-only.
84  libmesh_example_requires (3 == LIBMESH_DIM, "3D support");
85 
86  // Our input mesh here is in ExodusII format
87 #ifndef LIBMESH_HAVE_EXODUS_API
88  libmesh_example_requires (false, "ExodusII support");
89 #endif
90 
91 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
92  libmesh_example_requires (false, "second derivatives enabled");
93 #endif
94 
95  // This example does a bunch of linear algebra during assembly, and
96  // therefore requires Eigen.
97 #ifndef LIBMESH_HAVE_EIGEN
98  libmesh_example_requires(false, "--enable-eigen");
99 #endif
100 
101  // This example converts between ExodusII and XDR files, therefore
102  // it requires XDR support in libmesh.
103 #ifndef LIBMESH_HAVE_XDR
104  libmesh_example_requires (false, "XDR support");
105 #endif
106 
107  // Read the "distributed_load" flag from the command line
108  GetPot command_line (argc, argv);
109  int distributed_load = 0;
110  if (command_line.search(1, "-distributed_load"))
111  distributed_load = command_line.next(distributed_load);
112 
113  {
114  Mesh mesh (init.comm(), 3);
115 
116  // To confirm that both ExodusII and Xdr formats work for shell
117  // meshes, we read in cylinder.exo, then write out cylinder.xdr,
118  // then read in cylinder.exo again below and use that for the rest
119  // of the example.
120  mesh.read("cylinder.exo");
121  mesh.write("cylinder.xdr");
122  }
123  Mesh mesh (init.comm(), 3);
124  mesh.read("cylinder.xdr");
125 
126  // Print information about the mesh to the screen.
127  mesh.print_info();
128 
129  // Create an equation systems object.
130  EquationSystems equation_systems (mesh);
131 
132  // Declare the system and its variables.
133  // Create a linear implicit system named "Shell".
134  LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem> ("Shell");
135 
136  // Add the three displacement variables "u", "v", "w",
137  // and the three rotational variables "theta_x", "theta_y", "theta_z".
138  // All variables are second order.
139  system.add_variable ("u",SECOND,LAGRANGE);
140  system.add_variable ("v",SECOND,LAGRANGE);
141  system.add_variable ("w",SECOND,LAGRANGE);
142  system.add_variable ("theta_x",SECOND,LAGRANGE);
143  system.add_variable ("theta_y",SECOND,LAGRANGE);
144  system.add_variable ("theta_z",SECOND,LAGRANGE);
145 
146  // Give the system a pointer to the matrix and rhs assembly
147  // function.
149 
150  // Use the parameters of the equation systems object to
151  // tell the shell system about the material properties, the
152  // shell thickness, and the external load.
153  const Real h = 0.03;
154  const Real E = 3e10;
155  const Real nu = 0.3;
156  const Real q = 1;
157  equation_systems.parameters.set<Real> ("thickness") = h;
158  equation_systems.parameters.set<Real> ("young's modulus") = E;
159  equation_systems.parameters.set<Real> ("poisson ratio") = nu;
160  equation_systems.parameters.set<Real> ("point load") = q;
161  equation_systems.parameters.set<bool>("distributed load") = (distributed_load != 0);
162 
163  // Dirichlet conditions for the pinched cylinder problem.
164  // Only one 8th of the cylinder is considered using symmetry considerations.
165  // The cylinder longitudinal axis is the y-axis.
166  // The four corners of the surface are named A(3,0,0), B(3,3,0), C(0,3,3), D(0,0,3).
167  // The point load (pinch) is applied at C in the -z direction.
168  // Edge AD is the actual edge of the cylinder and is rigid in the xz-plane.
169  // Other edges have symmetric boundary conditions.
170 
171  // AB w, theta_x, theta_y
172  {
173  std::set<boundary_id_type> boundary_ids;
174  boundary_ids.insert(7);
175  unsigned int variables[] = {2, 3, 4};
176  ZeroFunction<> zf;
177 
178  // Most DirichletBoundary users will want to supply a "locally
179  // indexed" functor
180  DirichletBoundary dirichlet_bc
181  (boundary_ids,
182  std::vector<unsigned int>(variables, variables+3), zf,
184  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
185  }
186  // BC v, theta_x, theta_z
187  {
188  std::set<boundary_id_type> boundary_ids;
189  boundary_ids.insert(8);
190  unsigned int variables[] = {1, 3, 5};
191  ZeroFunction<> zf;
192 
193  DirichletBoundary dirichlet_bc
194  (boundary_ids,
195  std::vector<unsigned int>(variables, variables+3), zf,
197  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
198  }
199  // CD u, theta_y, theta_z
200  {
201  std::set<boundary_id_type> boundary_ids;
202  boundary_ids.insert(9);
203  unsigned int variables[] = {0, 4, 5};
204  ZeroFunction<> zf;
205 
206  DirichletBoundary dirichlet_bc
207  (boundary_ids,
208  std::vector<unsigned int>(variables, variables+3), zf,
210  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
211  }
212  // AD u, w, theta_y
213  {
214  std::set<boundary_id_type> boundary_ids;
215  boundary_ids.insert(10);
216  unsigned int variables[] = {0, 2, 4};
217  ZeroFunction<> zf;
218 
219  DirichletBoundary dirichlet_bc
220  (boundary_ids,
221  std::vector<unsigned int>(variables, variables+3), zf,
223  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
224  }
225 
226  // Initialize the data structures for the equation system.
227  equation_systems.init();
228 
229  // Print information about the system to the screen.
230  equation_systems.print_info();
231 
232  // This example can be run with EigenSparseLinearSolvers, but it
233  // only works with either the CG or SPARSELU types, and SparseLU
234  // turns out to be faster.
237 
238  // Solve the linear system.
239  system.solve();
240 
241  // After solving the system, write the solution to an
242  // ExodusII output file ready for import in, e.g.,
243  // Paraview.
244  ExodusII_IO(mesh).write_equation_systems ("out.e", equation_systems);
245 
246  // Compare with analytical solution for point load
247  if (distributed_load==0)
248  {
249  // Find the node nearest point C.
250  Node * node_C = libmesh_nullptr;
251  Point point_C(0, 3, 3);
252  {
253  Real nearest_dist_sq = std::numeric_limits<Real>::max();
254 
255  // Find the closest local node. On a DistributedMesh we may
256  // not even know about the existence of closer non-local
257  // nodes.
258  for (const auto & node : mesh.local_node_ptr_range())
259  {
260  const Real dist_sq = (*node - point_C).norm_sq();
261  if (dist_sq < nearest_dist_sq)
262  {
263  nearest_dist_sq = dist_sq;
264  node_C = node;
265  }
266  }
267 
268  // Check with other processors to see if any found a closer node
269  unsigned int minrank = 0;
270  system.comm().minloc(nearest_dist_sq, minrank);
271 
272  // Broadcast the ID of the closest node, so every processor can
273  // see for certain whether they have it or not.
274  dof_id_type nearest_node_id;
275  if (system.processor_id() == minrank)
276  nearest_node_id = node_C->id();
277  system.comm().broadcast(nearest_node_id, minrank);
278  node_C = mesh.query_node_ptr(nearest_node_id);
279  }
280 
281  // Evaluate the z-displacement "w" at the node nearest C.
282  Number w = 0;
283 
284  // If we know about the closest node, and if we also own the DoFs
285  // on that node, then we can evaluate the solution at that node.
286  if (node_C)
287  {
288  const unsigned int w_var = system.variable_number ("w");
289  dof_id_type w_dof = node_C->dof_number (system.number(), w_var, 0);
290  if (w_dof >= system.get_dof_map().first_dof() &&
291  w_dof < system.get_dof_map().end_dof())
292  w = system.current_solution(w_dof);
293  }
294  system.comm().sum(w);
295 
296 
297  Number w_C_bar = -E*h*w/q;
298  const Real w_C_bar_analytic = 164.24;
299 
300  // Print the finite element solution and the analytic
301  // prediction to the screen.
302  libMesh::out << "z-displacement of the point C: " << w_C_bar << std::endl;
303  libMesh::out << "Analytic solution: " << w_C_bar_analytic << std::endl;
304 
305  // Evaluate the y-displacement "v" at point D. This time we'll
306  // evaluate at the exact point, not just the closest node.
307  Point point_D(0, 0, 3);
308  const unsigned int v_var = system.variable_number ("v");
309  Number v = system.point_value(v_var, point_D);
310 
311  Number v_D_bar = E*h*v/q;
312  const Real v_D_bar_analytic = 4.114;
313 
314  // Print the finite element solution and the analytic
315  // prediction to the screen.
316  libMesh::out << "y-displacement of the point D: " << v_D_bar << std::endl;
317  libMesh::out << "Analytic solution: " << v_D_bar_analytic << std::endl;
318  }
319 
320  // All done.
321  return 0;
322 }
323 
324 
325 
326 // We now define the matrix and rhs vector assembly function
327 // for the shell system.
329  const std::string & system_name)
330 {
331  // This example requires Eigen to actually work, but we should still
332  // let it compile and throw a runtime error if you don't.
333 #ifdef LIBMESH_HAVE_EIGEN
334  // It is a good idea to make sure we are assembling
335  // the proper system.
336  libmesh_assert_equal_to (system_name, "Shell");
337 
338  // Get a constant reference to the mesh object.
339  const MeshBase & mesh = es.get_mesh();
340  const unsigned int dim = mesh.mesh_dimension();
341 
342  // Get a reference to the shell system object.
343  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> (system_name);
344 
345  // Get the shell parameters that we need during assembly.
346  const Real h = es.parameters.get<Real> ("thickness");
347  const Real E = es.parameters.get<Real> ("young's modulus");
348  const Real nu = es.parameters.get<Real> ("poisson ratio");
349  const Real q = es.parameters.get<Real> ("point load");
350  const bool distributed_load = es.parameters.get<bool> ("distributed load");
351 
352  // The membrane elastic matrix.
353  Eigen::Matrix3d Hm;
354  Hm <<
355  1., nu, 0.,
356  nu, 1., 0.,
357  0., 0., 0.5 * (1-nu);
358  Hm *= h * E/(1-nu*nu);
359 
360  // The bending elastic matrix.
361  Eigen::Matrix3d Hf;
362  Hf <<
363  1., nu, 0.,
364  nu, 1., 0.,
365  0., 0., 0.5 * (1-nu);
366  Hf *= h*h*h/12 * E/(1-nu*nu);
367 
368  // The shear elastic matrices.
369  Eigen::Matrix2d Hc0 = Eigen::Matrix2d::Identity();
370  Hc0 *= h * 5./6*E/(2*(1+nu));
371 
372  Eigen::Matrix2d Hc1 = Eigen::Matrix2d::Identity();
373  Hc1 *= h*h*h/12 * 5./6*E/(2*(1+nu));
374 
375  // Get the Finite Element type, this will be
376  // the same for all variables.
377  FEType fe_type = system.variable_type (0);
378 
379  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
380  QGauss qrule (dim, SECOND); // Reduced integration, 2x2 gauss points (instead of 3x3 for full integration)
381  fe->attach_quadrature_rule (&qrule);
382 
383  // The element Jacobian * quadrature weight at each integration point.
384  const std::vector<Real> & JxW = fe->get_JxW();
385 
386  // The element shape function and its derivatives evaluated at the
387  // quadrature points.
388  const std::vector<RealGradient> & dxyzdxi = fe->get_dxyzdxi();
389  const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();
390  const std::vector<RealGradient> & d2xyzdxi2 = fe->get_d2xyzdxi2();
391  const std::vector<RealGradient> & d2xyzdeta2 = fe->get_d2xyzdeta2();
392  const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();
393  const std::vector<std::vector<Real>> & dphidxi = fe->get_dphidxi();
394  const std::vector<std::vector<Real>> & dphideta = fe->get_dphideta();
395  const std::vector<std::vector<Real>> & phi = fe->get_phi();
396 
397  // A reference to the DofMap object for this system. The DofMap
398  // object handles the index translation from node and element numbers
399  // to degree of freedom numbers.
400  const DofMap & dof_map = system.get_dof_map();
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<Eigen::Matrix3d> Qnode;
446  for (unsigned int i=0; i<elem->n_nodes(); ++i)
447  {
448  Eigen::Vector3d a1;
449  a1 << dxyzdxi[i](0), dxyzdxi[i](1), dxyzdxi[i](2);
450  Eigen::Vector3d a2;
451  a2 << dxyzdeta[i](0), dxyzdeta[i](1), dxyzdeta[i](2);
452  Eigen::Vector3d 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  Eigen::Matrix3d 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  Eigen::Matrix3d 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  Eigen::Vector3d a1;
496  a1 << dxyzdxi[qp](0), dxyzdxi[qp](1), dxyzdxi[qp](2);
497  Eigen::Vector3d a2;
498  a2 << dxyzdeta[qp](0), dxyzdeta[qp](1), dxyzdeta[qp](2);
499  Eigen::Vector3d n;
500  n = a1.cross(a2);
501  n /= n.norm();
502  Eigen::Matrix3d 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  Eigen::Matrix3d 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  Eigen::Matrix3d 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  Eigen::Matrix2d C0;
533  C0 = F0it.block<3,2>(0,0).transpose()*Q.block<3,2>(0,0);
534 
535  // Normal derivatives in reference coordinates
536  Eigen::Vector3d d2Xdxi2(d2xyzdxi2[qp](0), d2xyzdxi2[qp](1), d2xyzdxi2[qp](2));
537  Eigen::Vector3d d2Xdeta2(d2xyzdeta2[qp](0), d2xyzdeta2[qp](1), d2xyzdeta2[qp](2));
538  Eigen::Vector3d d2Xdxideta(d2xyzdxideta[qp](0), d2xyzdxideta[qp](1), d2xyzdxideta[qp](2));
539 
540  Eigen::Matrix2d b;
541  b <<
542  n.dot(d2Xdxi2), n.dot(d2Xdxideta),
543  n.dot(d2Xdxideta), n.dot(d2Xdeta2);
544 
545  Eigen::Vector3d dndxi = -b(0,0)*F0it.col(0) - b(0,1)*F0it.col(1);
546  Eigen::Vector3d dndeta = -b(1,0)*F0it.col(0) - b(1,1)*F0it.col(1);
547 
548  Eigen::Matrix2d 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  Eigen::Matrix2d 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  Eigen::MatrixXd B0I(3, 5);
567  B0I = Eigen::MatrixXd::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  Eigen::Vector2d V1i(-Q.col(0).dot(Qnode[i].col(1)),
577  Q.col(0).dot(Qnode[i].col(0)));
578 
579  Eigen::Vector2d V2i(-Q.col(1).dot(Qnode[i].col(1)),
580  Q.col(1).dot(Qnode[i].col(0)));
581 
582  Eigen::MatrixXd B1I(3,5);
583  B1I = Eigen::MatrixXd::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  Eigen::MatrixXd B2I(3,5);
594  B2I = Eigen::MatrixXd::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  Eigen::MatrixXd Bc0I(2,5);
602  Bc0I = Eigen::MatrixXd::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  Eigen::MatrixXd Bc1I(2,5);
610  Bc1I = Eigen::MatrixXd::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  Eigen::Vector2d BdxiI(dphidxi[i][qp],dphideta[i][qp]);
616  Eigen::Vector2d 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  Eigen::MatrixXd B0J(3,5);
626  B0J = Eigen::MatrixXd::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  Eigen::Vector2d V1j(-Q.col(0).dot(Qnode[j].col(1)),
636  Q.col(0).dot(Qnode[j].col(0)));
637 
638  Eigen::Vector2d V2j(-Q.col(1).dot(Qnode[j].col(1)),
639  Q.col(1).dot(Qnode[j].col(0)));
640 
641  Eigen::MatrixXd B1J(3,5);
642  B1J = Eigen::MatrixXd::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  Eigen::MatrixXd B2J(3,5);
653  B2J = Eigen::MatrixXd::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  Eigen::MatrixXd Bc0J(2, 5);
661  Bc0J = Eigen::MatrixXd::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  Eigen::MatrixXd Bc1J(2, 5);
669  Bc1J = Eigen::MatrixXd::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  Eigen::Vector2d BdxiJ(dphidxi[j][qp], dphideta[j][qp]);
675  Eigen::Vector2d 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  Eigen::MatrixXd 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  Eigen::MatrixXd full_local_KIJ(6, 6);
695  full_local_KIJ = Eigen::MatrixXd::Zero(6, 6);
696  full_local_KIJ.block<5,5>(0,0)=local_KIJ;
697 
698  // Drilling dof stiffness contribution
699  full_local_KIJ(5,5) = Hf(0,0)*JxW[qp]*BdI.transpose()*BdJ;
700 
701  // Transform the stiffness matrix to global coordinates
702  Eigen::MatrixXd global_KIJ(6,6);
703  Eigen::MatrixXd TI(6,6);
704  TI = Eigen::MatrixXd::Identity(6,6);
705  TI.block<3,3>(3,3) = Qnode[i].transpose();
706  Eigen::MatrixXd TJ(6,6);
707  TJ = Eigen::MatrixXd::Identity(6,6);
708  TJ.block<3,3>(3,3) = Qnode[j].transpose();
709  global_KIJ = TI.transpose()*full_local_KIJ*TJ;
710 
711  // Insert the components of the coupling stiffness
712  // matrix KIJ into the corresponding directional
713  // submatrices.
714  for (unsigned int k=0;k<6;k++)
715  for (unsigned int l=0;l<6;l++)
716  Ke_var[k][l](i,j) += global_KIJ(k,l);
717  }
718  }
719 
720  } // end of the quadrature point qp-loop
721 
722  if (distributed_load)
723  {
724  // Loop on shell faces
725  for (unsigned int shellface=0; shellface<2; shellface++)
726  {
727  std::vector<boundary_id_type> bids;
728  mesh.get_boundary_info().shellface_boundary_ids(elem, shellface, bids);
729 
730  for (std::size_t k=0; k<bids.size(); k++)
731  if (bids[k]==11) // sideset id for surface load
732  for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
733  for (unsigned int i=0; i<n_var_dofs; ++i)
734  Fe_w(i) -= JxW[qp] * phi[i][qp];
735  }
736  }
737 
738  // The element matrix is now built for this element.
739  // Add it to the global matrix.
740 
741  dof_map.constrain_element_matrix_and_vector (Ke,Fe,dof_indices);
742 
743  system.matrix->add_matrix (Ke, dof_indices);
744  system.rhs->add_vector (Fe, dof_indices);
745 
746  }
747 
748  if (!distributed_load)
749  {
750  //Adding point load to the RHS
751 
752  //Pinch position
753  Point C(0, 3, 3);
754 
755  //Finish assembling rhs so we can set one value
756  system.rhs->close();
757 
759  const MeshBase::const_node_iterator node_end = mesh.nodes_end();
760 
761  for ( ; nodeit!=node_end; ++nodeit)
762  {
763  Node & node = **nodeit;
764  if ((node-C).norm() < 1e-3)
765  system.rhs->set(node.dof_number(0, 2, 0), -q/4);
766  }
767  }
768 
769 #else
770  // Avoid compiler warnings
771  libmesh_ignore(es);
772  libmesh_ignore(system_name);
773 #endif // LIBMESH_HAVE_EIGEN
774 }
Number point_value(unsigned int var, const Point &p, const bool insist_on_success=true) const
Definition: system.C:2011
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
int main(int argc, char **argv)
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
double abs(double a)
This is the EquationSystems class.
A Node is like a Point, but with more information.
Definition: node.h:52
virtual SimpleRange< node_iterator > local_node_ptr_range()=0
ConstFunction that simply returns 0.
Definition: zero_function.h:35
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
This class provides a specific system class.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
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]...
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
const T & get(const std::string &) const
Definition: parameters.h:430
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
long double max(long double a, double b)
void assemble_shell(EquationSystems &es, const std::string &system_name)
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
virtual node_iterator nodes_begin()=0
Iterate over all the nodes in the Mesh.
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:535
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
void reposition(const unsigned int ioff, const unsigned int n)
Changes the location of the subvector in the parent vector.
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:577
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.
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
const DofMap & get_dof_map() const
Definition: system.h:2030
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
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:1795
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 NULL/0 as necessary.
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:1806
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
virtual void write(const std::string &name)=0
void set_solver_type(const SolverType st)
Sets the type of solver to use.
void broadcast(T &data, const unsigned int root_id=0) const
Take a local value and broadcast it to all processors.
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_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 libmesh_ignore(const T &)
Definition: assembly.h:127
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual node_iterator nodes_end()=0
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
unsigned int number() const
Definition: system.h:2006
T & set(const std::string &)
Definition: parameters.h:469
const Parallel::Communicator & comm() const
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
The definition of the const_node_iterator struct.
Definition: mesh_base.h:1548
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
const T_sys & get_system(const std::string &name) const
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 id() const
Definition: dof_object.h:632
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
unsigned int n_points() const
Definition: quadrature.h:113
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64
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:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
void minloc(T &r, unsigned int &min_id) const
Take a local variable and replace it with the minimum of it&#39;s values on all processors, returning the minimum rank of a processor which originally held the minimum value.