libMesh
Functions
miscellaneous_ex9.C File Reference

Go to the source code of this file.

Functions

void assemble_poisson (EquationSystems &es, const ElementSideMap &lower_to_upper)
 Assemble the system matrix and rhs vector. More...
 
int main (int argc, char **argv)
 

Function Documentation

void assemble_poisson ( EquationSystems es,
const ElementSideMap lower_to_upper 
)

Assemble the system matrix and rhs vector.

Definition at line 206 of file miscellaneous_ex9.C.

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, libMesh::Parameters::get(), libMesh::MeshBase::get_boundary_info(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::BoundaryInfo::has_boundary_id(), libMesh::FEInterface::inverse_map(), libMesh::libmesh_assert(), libmesh_nullptr, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), side, and libMesh::TypeVector< T >::size().

Referenced by main().

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
unsigned int dim
This class provides a specific system class.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
unsigned short int side
Definition: xdr_io.C:49
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
const T & get(const std::string &) const
Definition: parameters.h:430
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
Order default_quadrature_order() const
Definition: fe_type.h:332
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
const T_sys & get_system(const std::string &name) const
int main ( int  argc,
char **  argv 
)

Definition at line 97 of file miscellaneous_ex9.C.

References libMesh::DofMap::add_coupling_functor(), libMesh::DofMap::add_dirichlet_boundary(), libMesh::System::add_variable(), libMesh::System::assemble_before_solve, assemble_poisson(), libMesh::LibMeshInit::comm(), libMesh::default_solver_package(), libMesh::FIRST, libMesh::System::get_dof_map(), AugmentSparsityOnInterface::get_lower_to_upper(), libMesh::TriangleWrapper::init(), libMesh::LAGRANGE, libMesh::LOCAL_VARIABLE_ORDER, mesh, libMesh::out, libMesh::PETSC_SOLVERS, libMesh::MeshBase::read(), libMesh::Real, libMesh::LinearImplicitSystem::solve(), libMesh::MeshRefinement::uniformly_refine(), and libMesh::MeshOutput< MT >::write_equation_systems().

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 }
This is the EquationSystems class.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
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.
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
void assemble_poisson(EquationSystems &es, const ElementSideMap &lower_to_upper)
Assemble the system matrix and rhs vector.
MeshBase & mesh
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
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
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
SolverPackage default_solver_package()
Definition: libmesh.C:995
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
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
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
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.