libMesh
miscellaneous_ex11.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 
19 
20 // <h1>Miscellaneous Example 11 - Using Loop Subdivision Shell Elements</h1>
21 // \author Roman Vetter
22 // \author Norbert Stoop
23 // \date 2014
24 //
25 // This example demonstrates how subdivision surface shell elements
26 // are used, and how boundary conditions can be applied to them. To
27 // keep it simple, we solve the static deflection of a clamped, square,
28 // linearly elastic Kirchhoff-Love plate subject to a uniform
29 // load distribution. Refer to Cirak et al., Int. J. Numer. Meth.
30 // Engng. 2000; 47: 2039-2072, for a detailed description of what's
31 // implemented. In fact, this example follows that paper very closely.
32 
33 
34 // C++ include files that we need
35 #include <iostream>
36 
37 // LibMesh includes
38 #include "libmesh/libmesh.h"
39 #include "libmesh/replicated_mesh.h"
40 #include "libmesh/mesh_refinement.h"
41 #include "libmesh/mesh_modification.h"
42 #include "libmesh/mesh_tools.h"
43 #include "libmesh/linear_implicit_system.h"
44 #include "libmesh/equation_systems.h"
45 #include "libmesh/fe.h"
46 #include "libmesh/quadrature.h"
47 #include "libmesh/node.h"
48 #include "libmesh/elem.h"
49 #include "libmesh/dof_map.h"
50 #include "libmesh/vector_value.h"
51 #include "libmesh/tensor_value.h"
52 #include "libmesh/dense_matrix.h"
53 #include "libmesh/dense_submatrix.h"
54 #include "libmesh/dense_vector.h"
55 #include "libmesh/dense_subvector.h"
56 #include "libmesh/sparse_matrix.h"
57 #include "libmesh/numeric_vector.h"
58 #include "libmesh/vtk_io.h"
59 #include "libmesh/exodusII_io.h"
60 
61 // These are the include files typically needed for subdivision elements.
62 #include "libmesh/face_tri3_subdivision.h"
63 #include "libmesh/mesh_subdivision_support.h"
64 
65 // Bring in everything from the libMesh namespace
66 using namespace libMesh;
67 
68 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
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 #endif
75 
76 // Begin the main program.
77 int main (int argc, char ** argv)
78 {
79  // Initialize libMesh.
80  LibMeshInit init (argc, argv);
81 
82  // This example requires a linear solver package.
83  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
84  "--enable-petsc, --enable-trilinos, or --enable-eigen");
85 
86  // Skip this 3D example if libMesh was compiled as 1D/2D-only.
87  libmesh_example_requires (3 == LIBMESH_DIM, "3D support");
88 
89  // Skip this example without --enable-node-valence
90 #ifndef LIBMESH_ENABLE_NODE_VALENCE
91  libmesh_example_requires (false, "--enable-node-valence");
92 #endif
93 
94  // Skip this example without --enable-amr; requires MeshRefinement
95 #ifndef LIBMESH_ENABLE_AMR
96  libmesh_example_requires(false, "--enable-amr");
97 #else
98 
99  // Skip this example without --enable-second; requires d2phi
100 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
101  libmesh_example_requires(false, "--enable-second");
102 #else
103 
104  // Create a 2D mesh distributed across the default MPI communicator.
105  // Subdivision surfaces do not appear to work with DistributedMesh yet.
106  ReplicatedMesh mesh (init.comm(), 2);
107 
108  // Read the coarse square mesh.
109  mesh.read ("square_mesh.off");
110 
111  // Resize the square plate to edge length L.
112  const Real L = 100.;
114 
115  // Quadrisect the mesh triangles a few times to obtain a
116  // finer mesh. Subdivision surface elements require the
117  // refinement data to be removed afterward.
118  MeshRefinement mesh_refinement (mesh);
119  mesh_refinement.uniformly_refine (3);
121 
122  // Write the mesh before the ghost elements are added.
123 #if defined(LIBMESH_HAVE_VTK)
124  VTKIO(mesh).write ("without_ghosts.pvtu");
125 #endif
126 #if defined(LIBMESH_HAVE_EXODUS_API)
127  ExodusII_IO(mesh).write ("without_ghosts.e");
128 #endif
129 
130  // Print information about the triangulated mesh to the screen.
131  mesh.print_info();
132 
133  // Turn the triangulated mesh into a subdivision mesh
134  // and add an additional row of "ghost" elements around
135  // it in order to complete the extended local support of
136  // the triangles at the boundaries. If the second
137  // argument is set to true, the outermost existing
138  // elements are converted into ghost elements, and the
139  // actual physical mesh is thus getting smaller.
141 
142  // Print information about the subdivision mesh to the screen.
143  mesh.print_info();
144 
145  // Write the mesh with the ghost elements added.
146  // Compare this to the original mesh to see the difference.
147 #if defined(LIBMESH_HAVE_VTK)
148  VTKIO(mesh).write ("with_ghosts.pvtu");
149 #endif
150 #if defined(LIBMESH_HAVE_EXODUS_API)
151  ExodusII_IO(mesh).write ("with_ghosts.e");
152 #endif
153 
154  // Create an equation systems object.
155  EquationSystems equation_systems (mesh);
156 
157  // Declare the system and its variables.
158  // Create a linear implicit system named "Shell".
159  LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem> ("Shell");
160 
161  // Add the three translational deformation variables
162  // "u", "v", "w" to "Shell". Since subdivision shell
163  // elements meet the C1-continuity requirement, no
164  // rotational or other auxiliary variables are needed.
165  // Loop Subdivision Elements are always interpolated
166  // by quartic box splines, hence the order must always
167  // be FOURTH.
168  system.add_variable ("u", FOURTH, SUBDIVISION);
169  system.add_variable ("v", FOURTH, SUBDIVISION);
170  system.add_variable ("w", FOURTH, SUBDIVISION);
171 
172  // Give the system a pointer to the matrix and rhs assembly
173  // function.
175 
176  // Use the parameters of the equation systems object to
177  // tell the shell system about the material properties, the
178  // shell thickness, and the external load.
179  const Real h = 1.;
180  const Real E = 1.e7;
181  const Real nu = 0.;
182  const Real q = 1.;
183  equation_systems.parameters.set<Real> ("thickness") = h;
184  equation_systems.parameters.set<Real> ("young's modulus") = E;
185  equation_systems.parameters.set<Real> ("poisson ratio") = nu;
186  equation_systems.parameters.set<Real> ("uniform load") = q;
187 
188  // Initialize the data structures for the equation system.
189  equation_systems.init();
190 
191  // Print information about the system to the screen.
192  equation_systems.print_info();
193 
194  // Solve the linear system.
195  system.solve();
196 
197  // After solving the system, write the solution to a VTK
198  // or ExodusII output file ready for import in, e.g.,
199  // Paraview.
200 #if defined(LIBMESH_HAVE_VTK)
201  VTKIO(mesh).write_equation_systems ("out.pvtu", equation_systems);
202 #endif
203 #if defined(LIBMESH_HAVE_EXODUS_API)
204  ExodusII_IO(mesh).write_equation_systems ("out.e", equation_systems);
205 #endif
206 
207  // Find the center node to measure the maximum deformation of the plate.
208  Node * center_node = 0;
209  Real nearest_dist_sq = mesh.point(0).norm_sq();
210  for (unsigned int nid=1; nid<mesh.n_nodes(); ++nid)
211  {
212  const Real dist_sq = mesh.point(nid).norm_sq();
213  if (dist_sq < nearest_dist_sq)
214  {
215  nearest_dist_sq = dist_sq;
216  center_node = mesh.node_ptr(nid);
217  }
218  }
219 
220  // Finally, we evaluate the z-displacement "w" at the center node.
221  const unsigned int w_var = system.variable_number ("w");
222  dof_id_type w_dof = center_node->dof_number (system.number(), w_var, 0);
223  Number w = 0;
224  if (w_dof >= system.get_dof_map().first_dof() &&
225  w_dof < system.get_dof_map().end_dof())
226  w = system.current_solution(w_dof);
227  system.comm().sum(w);
228 
229 
230  // The analytic solution for the maximum displacement of
231  // a clamped square plate in pure bending, from Taylor,
232  // Govindjee, Commun. Numer. Meth. Eng. 20, 757-765, 2004.
233  const Real D = E * h*h*h / (12*(1-nu*nu));
234  const Real w_analytic = 0.001265319 * L*L*L*L * q / D;
235 
236  // Print the finite element solution and the analytic
237  // prediction of the maximum displacement of the clamped
238  // square plate to the screen.
239  libMesh::out << "z-displacement of the center point: " << w << std::endl;
240  libMesh::out << "Analytic solution for pure bending: " << w_analytic << std::endl;
241 
242 #endif // #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
243 
244 #endif // #ifdef LIBMESH_ENABLE_AMR
245 
246  // All done.
247  return 0;
248 }
249 
250 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
251 
252 // We now define the matrix and rhs vector assembly function
253 // for the shell system. This function implements the
254 // linear Kirchhoff-Love theory for thin shells. At the
255 // end we also take into account the boundary conditions
256 // here, using the penalty method.
258  const std::string & libmesh_dbg_var(system_name))
259 {
260  // It is a good idea to make sure we are assembling
261  // the proper system.
262  libmesh_assert_equal_to (system_name, "Shell");
263 
264  // Get a constant reference to the mesh object.
265  const MeshBase & mesh = es.get_mesh();
266 
267  // Get a reference to the shell system object.
268  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> ("Shell");
269 
270  // Get the shell parameters that we need during assembly.
271  const Real h = es.parameters.get<Real> ("thickness");
272  const Real E = es.parameters.get<Real> ("young's modulus");
273  const Real nu = es.parameters.get<Real> ("poisson ratio");
274  const Real q = es.parameters.get<Real> ("uniform load");
275 
276  // Compute the membrane stiffness K and the bending
277  // rigidity D from these parameters.
278  const Real K = E * h / (1-nu*nu);
279  const Real D = E * h*h*h / (12*(1-nu*nu));
280 
281  // Numeric ids corresponding to each variable in the system.
282  const unsigned int u_var = system.variable_number ("u");
283  const unsigned int v_var = system.variable_number ("v");
284  const unsigned int w_var = system.variable_number ("w");
285 
286  // Get the Finite Element type for "u". Note this will be
287  // the same as the type for "v" and "w".
288  FEType fe_type = system.variable_type (u_var);
289 
290  // Build a Finite Element object of the specified type.
291  UniquePtr<FEBase> fe (FEBase::build(2, fe_type));
292 
293  // A Gauss quadrature rule for numerical integration.
294  // For subdivision shell elements, a single Gauss point per
295  // element is sufficient, hence we use extraorder = 0.
296  const int extraorder = 0;
297  UniquePtr<QBase> qrule (fe_type.default_quadrature_rule (2, extraorder));
298 
299  // Tell the finite element object to use our quadrature rule.
300  fe->attach_quadrature_rule (qrule.get());
301 
302  // The element Jacobian * quadrature weight at each integration point.
303  const std::vector<Real> & JxW = fe->get_JxW();
304 
305  // The surface tangents in both directions at the quadrature points.
306  const std::vector<RealGradient> & dxyzdxi = fe->get_dxyzdxi();
307  const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();
308 
309  // The second partial derivatives at the quadrature points.
310  const std::vector<RealGradient> & d2xyzdxi2 = fe->get_d2xyzdxi2();
311  const std::vector<RealGradient> & d2xyzdeta2 = fe->get_d2xyzdeta2();
312  const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();
313 
314  // The element shape function and its derivatives evaluated at the
315  // quadrature points.
316  const std::vector<std::vector<Real>> & phi = fe->get_phi();
317  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
318  const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
319 
320  // A reference to the DofMap object for this system. The DofMap
321  // object handles the index translation from node and element numbers
322  // to degree of freedom numbers.
323  const DofMap & dof_map = system.get_dof_map();
324 
325  // Define data structures to contain the element stiffness matrix
326  // and right-hand-side vector contribution. Following
327  // basic finite element terminology we will denote these
328  // "Ke" and "Fe".
331 
333  Kuu(Ke), Kuv(Ke), Kuw(Ke),
334  Kvu(Ke), Kvv(Ke), Kvw(Ke),
335  Kwu(Ke), Kwv(Ke), Kww(Ke);
336 
338  Fu(Fe),
339  Fv(Fe),
340  Fw(Fe);
341 
342  // This vector will hold the degree of freedom indices for
343  // the element. These define where in the global system
344  // the element degrees of freedom get mapped.
345  std::vector<dof_id_type> dof_indices;
346  std::vector<dof_id_type> dof_indices_u;
347  std::vector<dof_id_type> dof_indices_v;
348  std::vector<dof_id_type> dof_indices_w;
349 
350  // Now we will loop over all the elements in the mesh. We will
351  // compute the element matrix and right-hand-side contribution.
352  for (const auto & elem : mesh.active_local_element_ptr_range())
353  {
354  // The ghost elements at the boundaries need to be excluded
355  // here, as they don't belong to the physical shell,
356  // but serve for a proper boundary treatment only.
357  libmesh_assert_equal_to (elem->type(), TRI3SUBDIVISION);
358  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *> (elem);
359  if (sd_elem->is_ghost())
360  continue;
361 
362  // Get the degree of freedom indices for the
363  // current element. These define where in the global
364  // matrix and right-hand-side this element will
365  // contribute to.
366  dof_map.dof_indices (elem, dof_indices);
367  dof_map.dof_indices (elem, dof_indices_u, u_var);
368  dof_map.dof_indices (elem, dof_indices_v, v_var);
369  dof_map.dof_indices (elem, dof_indices_w, w_var);
370 
371  const std::size_t n_dofs = dof_indices.size();
372  const std::size_t n_u_dofs = dof_indices_u.size();
373  const std::size_t n_v_dofs = dof_indices_v.size();
374  const std::size_t n_w_dofs = dof_indices_w.size();
375 
376  // Compute the element-specific data for the current
377  // element. This involves computing the location of the
378  // quadrature points and the shape functions
379  // (phi, dphi, d2phi) for the current element.
380  fe->reinit (elem);
381 
382  // Zero the element matrix and right-hand side before
383  // summing them. We use the resize member here because
384  // the number of degrees of freedom might have changed from
385  // the last element.
386  Ke.resize (n_dofs, n_dofs);
387  Fe.resize (n_dofs);
388 
389  // Reposition the submatrices... The idea is this:
390  //
391  // - - - -
392  // | Kuu Kuv Kuw | | Fu |
393  // Ke = | Kvu Kvv Kvw |; Fe = | Fv |
394  // | Kwu Kwv Kww | | Fw |
395  // - - - -
396  //
397  // The DenseSubMatrix.reposition () member takes the
398  // (row_offset, column_offset, row_size, column_size).
399  //
400  // Similarly, the DenseSubVector.reposition () member
401  // takes the (row_offset, row_size)
402  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
403  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
404  Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
405 
406  Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
407  Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
408  Kvw.reposition (v_var*n_v_dofs, w_var*n_v_dofs, n_v_dofs, n_w_dofs);
409 
410  Kwu.reposition (w_var*n_w_dofs, u_var*n_w_dofs, n_w_dofs, n_u_dofs);
411  Kwv.reposition (w_var*n_w_dofs, v_var*n_w_dofs, n_w_dofs, n_v_dofs);
412  Kww.reposition (w_var*n_w_dofs, w_var*n_w_dofs, n_w_dofs, n_w_dofs);
413 
414  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
415  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
416  Fw.reposition (w_var*n_u_dofs, n_w_dofs);
417 
418  // Now we will build the element matrix and right-hand-side.
419  for (unsigned int qp=0; qp<qrule->n_points(); ++qp)
420  {
421  // First, we compute the external force resulting
422  // from a load q distributed uniformly across the plate.
423  // Since the load is supposed to be transverse to the plate,
424  // it affects the z-direction, i.e. the "w" variable.
425  for (unsigned int i=0; i<n_u_dofs; ++i)
426  Fw(i) += JxW[qp] * phi[i][qp] * q;
427 
428  // Next, we assemble the stiffness matrix. This is only valid
429  // for the linear theory, i.e., for small deformations, where
430  // reference and deformed surface metrics are indistinguishable.
431 
432  // Get the three surface basis vectors.
433  const RealVectorValue & a1 = dxyzdxi[qp];
434  const RealVectorValue & a2 = dxyzdeta[qp];
435  RealVectorValue a3 = a1.cross(a2);
436  const Real jac = a3.norm(); // the surface Jacobian
437  libmesh_assert_greater (jac, 0);
438  a3 /= jac; // the shell director a3 is normalized to unit length
439 
440  // Get the derivatives of the surface tangents.
441  const RealVectorValue & a11 = d2xyzdxi2[qp];
442  const RealVectorValue & a22 = d2xyzdeta2[qp];
443  const RealVectorValue & a12 = d2xyzdxideta[qp];
444 
445  // Compute the three covariant components of the first
446  // fundamental form of the surface.
447  const RealVectorValue a(a1*a1, a2*a2, a1*a2);
448 
449  // The elastic H matrix in Voigt's notation, computed from the
450  // covariant components of the first fundamental form rather
451  // than the contravariant components, exploiting that the
452  // contravariant first fundamental form is the inverse of the
453  // covariant first fundamental form (hence the determinant etc.).
454  RealTensorValue H;
455  H(0,0) = a(1) * a(1);
456  H(0,1) = H(1,0) = nu * a(1) * a(0) + (1-nu) * a(2) * a(2);
457  H(0,2) = H(2,0) = -a(1) * a(2);
458  H(1,1) = a(0) * a(0);
459  H(1,2) = H(2,1) = -a(0) * a(2);
460  H(2,2) = 0.5 * ((1-nu) * a(1) * a(0) + (1+nu) * a(2) * a(2));
461  const Real det = a(0) * a(1) - a(2) * a(2);
462  libmesh_assert_not_equal_to (det * det, 0);
463  H /= det * det;
464 
465  // Precompute come cross products for the bending part below.
466  const RealVectorValue a11xa2 = a11.cross(a2);
467  const RealVectorValue a22xa2 = a22.cross(a2);
468  const RealVectorValue a12xa2 = a12.cross(a2);
469  const RealVectorValue a1xa11 = a1.cross(a11);
470  const RealVectorValue a1xa22 = a1.cross(a22);
471  const RealVectorValue a1xa12 = a1.cross(a12);
472  const RealVectorValue a2xa3 = a2.cross(a3);
473  const RealVectorValue a3xa1 = a3.cross(a1);
474 
475  // Loop over all pairs of nodes I,J.
476  for (unsigned int i=0; i<n_u_dofs; ++i)
477  {
478  for (unsigned int j=0; j<n_u_dofs; ++j)
479  {
480  // The membrane strain matrices in Voigt's notation.
481  RealTensorValue MI, MJ;
482  for (unsigned int k=0; k<3; ++k)
483  {
484  MI(0,k) = dphi[i][qp](0) * a1(k);
485  MI(1,k) = dphi[i][qp](1) * a2(k);
486  MI(2,k) = dphi[i][qp](1) * a1(k)
487  + dphi[i][qp](0) * a2(k);
488 
489  MJ(0,k) = dphi[j][qp](0) * a1(k);
490  MJ(1,k) = dphi[j][qp](1) * a2(k);
491  MJ(2,k) = dphi[j][qp](1) * a1(k)
492  + dphi[j][qp](0) * a2(k);
493  }
494 
495  // The bending strain matrices in Voigt's notation.
496  RealTensorValue BI, BJ;
497  for (unsigned int k=0; k<3; ++k)
498  {
499  const Real term_ik = dphi[i][qp](0) * a2xa3(k)
500  + dphi[i][qp](1) * a3xa1(k);
501  BI(0,k) = -d2phi[i][qp](0,0) * a3(k)
502  +(dphi[i][qp](0) * a11xa2(k)
503  + dphi[i][qp](1) * a1xa11(k)
504  + (a3*a11) * term_ik) / jac;
505  BI(1,k) = -d2phi[i][qp](1,1) * a3(k)
506  +(dphi[i][qp](0) * a22xa2(k)
507  + dphi[i][qp](1) * a1xa22(k)
508  + (a3*a22) * term_ik) / jac;
509  BI(2,k) = 2 * (-d2phi[i][qp](0,1) * a3(k)
510  +(dphi[i][qp](0) * a12xa2(k)
511  + dphi[i][qp](1) * a1xa12(k)
512  + (a3*a12) * term_ik) / jac);
513 
514  const Real term_jk = dphi[j][qp](0) * a2xa3(k)
515  + dphi[j][qp](1) * a3xa1(k);
516  BJ(0,k) = -d2phi[j][qp](0,0) * a3(k)
517  +(dphi[j][qp](0) * a11xa2(k)
518  + dphi[j][qp](1) * a1xa11(k)
519  + (a3*a11) * term_jk) / jac;
520  BJ(1,k) = -d2phi[j][qp](1,1) * a3(k)
521  +(dphi[j][qp](0) * a22xa2(k)
522  + dphi[j][qp](1) * a1xa22(k)
523  + (a3*a22) * term_jk) / jac;
524  BJ(2,k) = 2 * (-d2phi[j][qp](0,1) * a3(k)
525  +(dphi[j][qp](0) * a12xa2(k)
526  + dphi[j][qp](1) * a1xa12(k)
527  + (a3*a12) * term_jk) / jac);
528  }
529 
530  // The total stiffness matrix coupling the nodes
531  // I and J is a sum of membrane and bending
532  // contributions according to the following formula.
533  const RealTensorValue KIJ = JxW[qp] * K * MI.transpose() * H * MJ
534  + JxW[qp] * D * BI.transpose() * H * BJ;
535 
536  // Insert the components of the coupling stiffness
537  // matrix KIJ into the corresponding directional
538  // submatrices.
539  Kuu(i,j) += KIJ(0,0);
540  Kuv(i,j) += KIJ(0,1);
541  Kuw(i,j) += KIJ(0,2);
542 
543  Kvu(i,j) += KIJ(1,0);
544  Kvv(i,j) += KIJ(1,1);
545  Kvw(i,j) += KIJ(1,2);
546 
547  Kwu(i,j) += KIJ(2,0);
548  Kwv(i,j) += KIJ(2,1);
549  Kww(i,j) += KIJ(2,2);
550  }
551  }
552 
553  } // end of the quadrature point qp-loop
554 
555  // The element matrix and right-hand-side are now built
556  // for this element. Add them to the global matrix and
557  // right-hand-side vector. The NumericMatrix::add_matrix()
558  // and NumericVector::add_vector() members do this for us.
559  system.matrix->add_matrix (Ke, dof_indices);
560  system.rhs->add_vector (Fe, dof_indices);
561  } // end of non-ghost element loop
562 
563  // Next, we apply the boundary conditions. In this case,
564  // all boundaries are clamped by the penalty method, using
565  // the special "ghost" nodes along the boundaries. Note
566  // that there are better ways to implement boundary conditions
567  // for subdivision shells. We use the simplest way here,
568  // which is known to be overly restrictive and will lead to
569  // a slightly too small deformation of the plate.
570  for (const auto & elem : mesh.active_local_element_ptr_range())
571  {
572  // For the boundary conditions, we only need to loop over
573  // the ghost elements.
574  libmesh_assert_equal_to (elem->type(), TRI3SUBDIVISION);
575  const Tri3Subdivision * gh_elem = static_cast<const Tri3Subdivision *> (elem);
576  if (!gh_elem->is_ghost())
577  continue;
578 
579  // Find the side which is part of the physical plate boundary,
580  // that is, the boundary of the original mesh without ghosts.
581  for (auto s : elem->side_index_range())
582  {
583  const Tri3Subdivision * nb_elem = static_cast<const Tri3Subdivision *> (elem->neighbor_ptr(s));
584  if (nb_elem == libmesh_nullptr || nb_elem->is_ghost())
585  continue;
586 
587  /*
588  * Determine the four nodes involved in the boundary
589  * condition treatment of this side. The MeshTools::Subdiv
590  * namespace provides lookup tables next and prev
591  * for an efficient determination of the next and previous
592  * nodes of an element, respectively.
593  *
594  * n4
595  * / \
596  * / gh \
597  * n2 ---- n3
598  * \ nb /
599  * \ /
600  * n1
601  */
602  const Node * nodes [4]; // n1, n2, n3, n4
603  nodes[1] = gh_elem->node_ptr(s); // n2
604  nodes[2] = gh_elem->node_ptr(MeshTools::Subdivision::next[s]); // n3
605  nodes[3] = gh_elem->node_ptr(MeshTools::Subdivision::prev[s]); // n4
606 
607  // The node in the interior of the domain, n1, is the
608  // hardest to find. Walk along the edges of element nb until
609  // we have identified it.
610  unsigned int n_int = 0;
611  nodes[0] = nb_elem->node_ptr(0);
612  while (nodes[0]->id() == nodes[1]->id() || nodes[0]->id() == nodes[2]->id())
613  nodes[0] = nb_elem->node_ptr(++n_int);
614 
615  // The penalty value. \f$ \frac{1}{\epsilon} \f$
616  const Real penalty = 1.e10;
617 
618  // With this simple method, clamped boundary conditions are
619  // obtained by penalizing the displacements of all four nodes.
620  // This ensures that the displacement field vanishes on the
621  // boundary side s.
622  for (unsigned int n=0; n<4; ++n)
623  {
624  const dof_id_type u_dof = nodes[n]->dof_number (system.number(), u_var, 0);
625  const dof_id_type v_dof = nodes[n]->dof_number (system.number(), v_var, 0);
626  const dof_id_type w_dof = nodes[n]->dof_number (system.number(), w_var, 0);
627  system.matrix->add (u_dof, u_dof, penalty);
628  system.matrix->add (v_dof, v_dof, penalty);
629  system.matrix->add (w_dof, w_dof, penalty);
630  }
631  }
632  } // end of ghost element loop
633 }
634 
635 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
A Node is like a Point, but with more information.
Definition: node.h:52
virtual const Point & point(const dof_id_type i) const =0
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
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
static const unsigned int prev[3]
A lookup table for the decrement modulo 3 operation, for iterating through the three nodes per elemen...
Real norm() const
Definition: type_vector.h:909
MeshBase & mesh
void assemble_shell(EquationSystems &es, const std::string &system_name)
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
virtual const Node * node_ptr(const dof_id_type i) const =0
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.
This class implements reading and writing meshes in the VTK format.
Definition: vtk_io.h:59
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
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.
TypeTensor< T > transpose() const
Definition: type_tensor.h:992
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.
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
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.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
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.
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition: type_vector.h:879
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
Real norm_sq() const
Definition: type_vector.h:940
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
virtual void write(const std::string &) libmesh_override
Output the mesh without solutions to a .pvtu file.
static const unsigned int next[3]
A lookup table for the increment modulo 3 operation, for iterating through the three nodes per elemen...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
unsigned int number() const
Definition: system.h:2006
virtual void write(const std::string &fname) libmesh_override
This method implements writing a mesh to a specified file.
Definition: exodusII_io.C:789
T & set(const std::string &)
Definition: parameters.h:469
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
const Parallel::Communicator & comm() const
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
void prepare_subdivision_mesh(MeshBase &mesh, bool ghosted=false)
Prepares the mesh for use with subdivision elements.
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
void reposition(const unsigned int ioff, const unsigned int joff, const unsigned int new_m, const unsigned int new_n)
Changes the location of the submatrix in the parent matrix.
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
virtual dof_id_type n_nodes() const =0
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.
int main(int argc, char **argv)
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
void flatten(MeshBase &mesh)
Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elem...
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:64
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
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.