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