libMesh
miscellaneous_ex9.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 9 - Implement an interface term to model
21 // a thermal "film resistance"</h1>
22 // \author David Knezevic
23 // \date 2013
24 //
25 // In this example we solve a Poisson problem, -\Laplacian u = f, with
26 // a non-standard interface condition on the domain interior which
27 // models a thermal "film resistance". The interface condition
28 // requires continuity of flux, and a jump in temperature proportional
29 // to the flux:
30 // \nabla u_1 \cdot n = \nabla u_2 \cdot n,
31 // u_1 - u_2 = R * \nabla u \cdot n
32 //
33 // To implement this PDE, we use two mesh subdomains, \Omega_1 and
34 // \Omega_2, with coincident boundaries, but which are not connected
35 // in the FE sense. Let \Gamma denote the coincident boundary. The
36 // term on \Gamma takes the form:
37 //
38 // 1/R * \int_\Gamma (u_1 - u_2) (v_1 - v_2) ds,
39 //
40 // where u_1, u_2 (resp. v_1, v_2) are the trial (resp. test)
41 // functions on either side of \Gamma. We implement this condition
42 // using C0 basis functions, but the "crack" in the mesh at \Gamma
43 // permits a discontinuity in the solution. We also impose a heat flux
44 // on the bottom surface of the mesh, and a zero Dirichlet condition
45 // on the top surface.
46 //
47 // In order to implement the interface condition, we need to augment
48 // the matrix sparsity pattern, which is handled by the class
49 // AugmentSparsityPatternOnInterface.
50 
51 
52 // C++ include files that we need
53 #include <iostream>
54 #include <limits>
55 
56 // libMesh includes
57 #include "libmesh/libmesh.h"
58 #include "libmesh/mesh.h"
59 #include "libmesh/mesh_generation.h"
60 #include "libmesh/mesh_refinement.h"
61 #include "libmesh/exodusII_io.h"
62 #include "libmesh/equation_systems.h"
63 #include "libmesh/fe.h"
64 #include "libmesh/quadrature_gauss.h"
65 #include "libmesh/dof_map.h"
66 #include "libmesh/sparse_matrix.h"
67 #include "libmesh/numeric_vector.h"
68 #include "libmesh/dense_matrix.h"
69 #include "libmesh/dense_vector.h"
70 #include "libmesh/getpot.h"
71 #include "libmesh/elem.h"
72 #include "libmesh/fe_interface.h"
73 #include "libmesh/boundary_info.h"
74 #include "libmesh/linear_implicit_system.h"
75 #include "libmesh/zero_function.h"
76 #include "libmesh/dirichlet_boundaries.h"
77 
78 // local includes
80 
81 // define the boundary IDs in the mesh
82 #define MIN_Z_BOUNDARY 1
83 #define MAX_Z_BOUNDARY 2
84 #define CRACK_BOUNDARY_LOWER 3
85 #define CRACK_BOUNDARY_UPPER 4
86 
87 // Bring in everything from the libMesh namespace
88 using namespace libMesh;
89 
94  const ElementSideMap & lower_to_upper);
95 
96 // The main program.
97 int main (int argc, char ** argv)
98 {
99  // Initialize libMesh.
100  LibMeshInit init (argc, argv);
101 
102  // This example uses an ExodusII input file
103 #ifndef LIBMESH_HAVE_EXODUS_API
104  libmesh_example_requires(false, "--enable-exodus");
105 #endif
106 
107  // The sparsity augmentation code requires PETSc
108  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
109 
110  // Skip this 3D example if libMesh was compiled as 1D or 2D-only.
111  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");
112 
113  GetPot command_line (argc, argv);
114 
115  Real R = 2.;
116  if (command_line.search(1, "-R"))
117  R = command_line.next(R);
118 
119  Mesh mesh(init.comm());
120 
121  EquationSystems equation_systems (mesh);
122 
123  LinearImplicitSystem & system =
124  equation_systems.add_system<LinearImplicitSystem> ("Poisson");
125  system.add_variable("u", FIRST, LAGRANGE);
126 
127  // We want to call assemble_poisson "manually" so that we can pass in
128  // lower_to_upper, hence set assemble_before_solve = false
129  system.assemble_before_solve = false;
130 
131  // Impose zero Dirichlet boundary condition on MAX_Z_BOUNDARY
132  std::set<boundary_id_type> boundary_ids;
133  boundary_ids.insert(MAX_Z_BOUNDARY);
134  std::vector<unsigned int> variables;
135  variables.push_back(0);
136  ZeroFunction<> zf;
137 
138  // Most DirichletBoundary users will want to supply a "locally
139  // indexed" functor
140  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
142  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
143 
144  // Attach an object to the DofMap that will augment the sparsity pattern
145  // due to the degrees-of-freedom on the "crack"
146  //
147  // By attaching this object *before* reading our mesh, we also
148  // ensure that the connected elements will not be deleted on a
149  // distributed mesh.
150  AugmentSparsityOnInterface augment_sparsity
151  (mesh, CRACK_BOUNDARY_LOWER, CRACK_BOUNDARY_UPPER);
152  system.get_dof_map().add_coupling_functor(augment_sparsity);
153 
154  mesh.read("miscellaneous_ex9.exo");
155 
156  equation_systems.init();
157  equation_systems.print_info();
158 
159  // Set the jump term coefficient, it will be used in assemble_poisson
160  equation_systems.parameters.set<Real>("R") = R;
161 
162  // Assemble and then solve
163  assemble_poisson(equation_systems,
164  augment_sparsity.get_lower_to_upper());
165  system.solve();
166 
167 #ifdef LIBMESH_HAVE_EXODUS_API
168  // Plot the solution
169  ExodusII_IO (mesh).write_equation_systems ("solution.exo",
170  equation_systems);
171 #endif
172 
173 #ifdef LIBMESH_ENABLE_AMR
174  // Possibly solve on a refined mesh next.
175  MeshRefinement mesh_refinement (mesh);
176  unsigned int n_refinements = 0;
177  if (command_line.search("-n_refinements"))
178  n_refinements = command_line.next(0);
179 
180  for (unsigned int r = 0; r != n_refinements; ++r)
181  {
182  std::cout << "Refining the mesh" << std::endl;
183 
184  mesh_refinement.uniformly_refine ();
185  equation_systems.reinit();
186 
187  assemble_poisson(equation_systems,
188  augment_sparsity.get_lower_to_upper());
189  system.solve();
190 
191 #ifdef LIBMESH_HAVE_EXODUS_API
192  // Plot the refined solution
193  std::ostringstream out;
194  out << "solution_" << r << ".exo";
196  equation_systems);
197 #endif
198 
199  }
200 
201 #endif
202 
203  return 0;
204 }
205 
207  const ElementSideMap & lower_to_upper)
208 {
209  const MeshBase & mesh = es.get_mesh();
210  const unsigned int dim = mesh.mesh_dimension();
211 
212  Real R = es.parameters.get<Real>("R");
213 
214  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Poisson");
215 
216  const DofMap & dof_map = system.get_dof_map();
217 
218  FEType fe_type = dof_map.variable_type(0);
219 
220  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
221  UniquePtr<FEBase> fe_elem_face (FEBase::build(dim, fe_type));
222  UniquePtr<FEBase> fe_neighbor_face (FEBase::build(dim, fe_type));
223 
224  QGauss qrule (dim, fe_type.default_quadrature_order());
225  QGauss qface(dim-1, fe_type.default_quadrature_order());
226 
227  fe->attach_quadrature_rule (&qrule);
228  fe_elem_face->attach_quadrature_rule (&qface);
229  fe_neighbor_face->attach_quadrature_rule (&qface);
230 
231  const std::vector<Real> & JxW = fe->get_JxW();
232  const std::vector<std::vector<Real>> & phi = fe->get_phi();
233  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
234 
235  const std::vector<Real> & JxW_face = fe_elem_face->get_JxW();
236 
237  const std::vector<Point> & qface_points = fe_elem_face->get_xyz();
238 
239  const std::vector<std::vector<Real>> & phi_face = fe_elem_face->get_phi();
240  const std::vector<std::vector<Real>> & phi_neighbor_face = fe_neighbor_face->get_phi();
241 
244 
249 
250  std::vector<dof_id_type> dof_indices;
251 
252  for (const auto & elem : mesh.active_local_element_ptr_range())
253  {
254  dof_map.dof_indices (elem, dof_indices);
255  const unsigned int n_dofs = dof_indices.size();
256 
257  fe->reinit (elem);
258 
259  Ke.resize (n_dofs, n_dofs);
260  Fe.resize (n_dofs);
261 
262  // Assemble element interior terms for the matrix
263  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
264  for (unsigned int i=0; i<n_dofs; i++)
265  for (unsigned int j=0; j<n_dofs; j++)
266  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
267 
268  // Boundary flux provides forcing in this example
269  {
270  for (auto side : elem->side_index_range())
271  if (elem->neighbor_ptr(side) == libmesh_nullptr)
272  {
273  if (mesh.get_boundary_info().has_boundary_id (elem, side, MIN_Z_BOUNDARY))
274  {
275  fe_elem_face->reinit(elem, side);
276 
277  for (unsigned int qp=0; qp<qface.n_points(); qp++)
278  for (std::size_t i=0; i<phi.size(); i++)
279  Fe(i) += JxW_face[qp] * phi_face[i][qp];
280  }
281 
282  }
283  }
284 
285  // Add boundary terms on the crack
286  {
287  for (auto side : elem->side_index_range())
288  if (elem->neighbor_ptr(side) == libmesh_nullptr)
289  {
290  // Found the lower side of the crack. Assemble terms due to lower and upper in here.
291  if (mesh.get_boundary_info().has_boundary_id (elem, side, CRACK_BOUNDARY_LOWER))
292  {
293  fe_elem_face->reinit(elem, side);
294 
295  ElementSideMap::const_iterator ltu_it =
296  lower_to_upper.find(std::make_pair(elem, side));
297  libmesh_assert(ltu_it != lower_to_upper.end());
298 
299  const Elem * neighbor = ltu_it->second;
300 
301  std::vector<Point> qface_neighbor_points;
302  FEInterface::inverse_map (elem->dim(), fe->get_fe_type(),
303  neighbor, qface_points, qface_neighbor_points);
304  fe_neighbor_face->reinit(neighbor, &qface_neighbor_points);
305 
306  std::vector<dof_id_type> neighbor_dof_indices;
307  dof_map.dof_indices (neighbor, neighbor_dof_indices);
308  const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
309 
310  Kne.resize (n_neighbor_dofs, n_dofs);
311  Ken.resize (n_dofs, n_neighbor_dofs);
312  Kee.resize (n_dofs, n_dofs);
313  Knn.resize (n_neighbor_dofs, n_neighbor_dofs);
314 
315  // Lower-to-lower coupling term
316  for (unsigned int qp=0; qp<qface.n_points(); qp++)
317  for (unsigned int i=0; i<n_dofs; i++)
318  for (unsigned int j=0; j<n_dofs; j++)
319  Kee(i,j) -= JxW_face[qp] * (1./R)*(phi_face[i][qp] * phi_face[j][qp]);
320 
321  // Lower-to-upper coupling term
322  for (unsigned int qp=0; qp<qface.n_points(); qp++)
323  for (unsigned int i=0; i<n_dofs; i++)
324  for (unsigned int j=0; j<n_neighbor_dofs; j++)
325  Ken(i,j) += JxW_face[qp] * (1./R)*(phi_face[i][qp] * phi_neighbor_face[j][qp]);
326 
327  // Upper-to-upper coupling term
328  for (unsigned int qp=0; qp<qface.n_points(); qp++)
329  for (unsigned int i=0; i<n_neighbor_dofs; i++)
330  for (unsigned int j=0; j<n_neighbor_dofs; j++)
331  Knn(i,j) -= JxW_face[qp] * (1./R)*(phi_neighbor_face[i][qp] * phi_neighbor_face[j][qp]);
332 
333  // Upper-to-lower coupling term
334  for (unsigned int qp=0; qp<qface.n_points(); qp++)
335  for (unsigned int i=0; i<n_neighbor_dofs; i++)
336  for (unsigned int j=0; j<n_dofs; j++)
337  Kne(i,j) += JxW_face[qp] * (1./R)*(phi_neighbor_face[i][qp] * phi_face[j][qp]);
338 
339  system.matrix->add_matrix(Kne, neighbor_dof_indices, dof_indices);
340  system.matrix->add_matrix(Ken, dof_indices, neighbor_dof_indices);
341  system.matrix->add_matrix(Kee, dof_indices);
342  system.matrix->add_matrix(Knn, neighbor_dof_indices);
343  }
344  }
345  }
346 
347  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
348 
349  system.matrix->add_matrix (Ke, dof_indices);
350  system.rhs->add_vector (Fe, dof_indices);
351  }
352 }
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
This is the EquationSystems class.
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
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
unsigned short int side
Definition: xdr_io.C:49
void assemble_poisson(EquationSystems &es, const ElementSideMap &lower_to_upper)
Assemble the system matrix and rhs vector.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
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
void add_coupling_functor(GhostingFunctor &coupling_functor)
Adds a functor which can specify coupling requirements for creation of sparse matrices.
Definition: dof_map.C:1794
const ElementSideMap & get_lower_to_upper() const
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 is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
int main(int argc, char **argv)
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
Order default_quadrature_order() const
Definition: fe_type.h:332
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:569
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
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
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
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
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
Real size() const
Definition: type_vector.h:898
bool assemble_before_solve
Flag which tells the system to whether or not to call the user assembly function during each call to ...
Definition: system.h:1477
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.
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
std::map< std::pair< const Elem *, unsigned char >, const Elem * > ElementSideMap
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.