libMesh
amr.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 #include "libmesh/coupling_matrix.h"
19 #include "libmesh/dense_matrix.h"
20 #include "libmesh/dense_vector.h"
21 #include "libmesh/dof_map.h"
22 #include "libmesh/elem.h"
23 #include "libmesh/equation_systems.h"
24 #include "libmesh/fe.h"
25 #include "libmesh/gmv_io.h"
26 #include "libmesh/libmesh.h"
27 #include "libmesh/linear_implicit_system.h"
28 #include "libmesh/mesh.h"
29 #include "libmesh/mesh_refinement.h"
30 #include "libmesh/numeric_vector.h"
31 #include "libmesh/quadrature_gauss.h"
32 #include "libmesh/sparse_matrix.h"
33 
34 
35 using namespace libMesh;
36 
37 
38 void assemble(EquationSystems & es,
39  const std::string & system_name);
40 
41 
42 
43 
44 
45 #ifdef LIBMESH_ENABLE_AMR
46 int main (int argc, char ** argv)
47 {
48  LibMeshInit init(argc, argv);
49 
50  if (argc < 4)
51  libmesh_error_msg("Usage: ./prog -d DIM filename");
52 
53  // Variables to get us started
54  const unsigned int dim = atoi(argv[2]);
55 
56  std::string meshname (argv[3]);
57 
58  // declare a mesh...
59  Mesh mesh(init.comm(), dim);
60 
61  // Read a mesh
62  mesh.read(meshname);
63 
64  GMVIO(mesh).write ("out_0.gmv");
65 
67 
68  MeshRefinement mesh_refinement (mesh);
69 
70  mesh_refinement.refine_and_coarsen_elements ();
71  mesh_refinement.uniformly_refine (2);
72 
73  mesh.print_info();
74 
75 
76  // Set up the equation system(s)
77  EquationSystems es (mesh);
78 
79  LinearImplicitSystem & primary =
80  es.add_system<LinearImplicitSystem>("primary");
81 
82  primary.add_variable ("U", FIRST);
83  primary.add_variable ("V", FIRST);
84 
85  primary.get_dof_map()._dof_coupling->resize(2);
86  (*primary.get_dof_map()._dof_coupling)(0,0) = 1;
87  (*primary.get_dof_map()._dof_coupling)(1,1) = 1;
88 
90 
91  es.init ();
92 
93  es.print_info ();
95 
96  // call the solver.
97  primary.solve ();
98 
99  GMVIO(mesh).write_equation_systems ("out_1.gmv",
100  es);
101 
102 
103 
104  // Refine uniformly
105  mesh_refinement.uniformly_refine (1);
106  es.reinit ();
107 
108  // Write out the projected solution
109  GMVIO(mesh).write_equation_systems ("out_2.gmv",
110  es);
111 
112  // Solve again. Output the refined solution
113  primary.solve ();
114  GMVIO(mesh).write_equation_systems ("out_3.gmv",
115  es);
116 
117  return 0;
118 }
119 #else
120 int main (int, char **)
121 {
122  return 1;
123 }
124 #endif // ENABLE_AMR
125 
126 
127 
128 
129 
130 
132  const std::string & libmesh_dbg_var(system_name))
133 {
134  libmesh_assert_equal_to (system_name, "primary");
135 
136  const MeshBase & mesh = es.get_mesh();
137  const unsigned int dim = mesh.mesh_dimension();
138 
139  // Also use a 3x3x3 quadrature rule (3D). Then tell the FE
140  // about the geometry of the problem and the quadrature rule
141  FEType fe_type (FIRST);
142 
143  UniquePtr<FEBase> fe(FEBase::build(dim, fe_type));
144  QGauss qrule(dim, FIFTH);
145 
146  fe->attach_quadrature_rule (&qrule);
147 
148  UniquePtr<FEBase> fe_face(FEBase::build(dim, fe_type));
149  QGauss qface(dim-1, FIFTH);
150 
151  fe_face->attach_quadrature_rule(&qface);
152 
153  LinearImplicitSystem & system =
154  es.get_system<LinearImplicitSystem>("primary");
155 
156 
157  // These are references to cell-specific data
158  const std::vector<Real> & JxW_face = fe_face->get_JxW();
159  const std::vector<Real> & JxW = fe->get_JxW();
160  const std::vector<Point> & q_point = fe->get_xyz();
161  const std::vector<std::vector<Real>> & phi = fe->get_phi();
162  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
163 
164  std::vector<dof_id_type> dof_indices_U;
165  std::vector<dof_id_type> dof_indices_V;
166  const DofMap & dof_map = system.get_dof_map();
167 
172 
173  Real vol=0., area=0.;
174 
175  for (const auto & elem : mesh.active_local_element_ptr_range())
176  {
177  // recompute the element-specific data for the current element
178  fe->reinit (elem);
179 
180 
181  //fe->print_info();
182 
183  dof_map.dof_indices(elem, dof_indices_U, 0);
184  dof_map.dof_indices(elem, dof_indices_V, 1);
185 
186  // zero the element matrix and vector
187  Kuu.resize (phi.size(),
188  phi.size());
189 
190  Kvv.resize (phi.size(),
191  phi.size());
192 
193  Fu.resize (phi.size());
194  Fv.resize (phi.size());
195 
196  // standard stuff... like in code 1.
197  for (unsigned int gp=0; gp<qrule.n_points(); gp++)
198  {
199  for (std::size_t i=0; i<phi.size(); ++i)
200  {
201  // this is tricky. ig is the _global_ dof index corresponding
202  // to the _global_ vertex number elem->node_id(i). Note that
203  // in general these numbers will not be the same (except for
204  // the case of one unknown per node on one subdomain) so
205  // we need to go through the dof_map
206 
207  const Real f = q_point[gp]*q_point[gp];
208  // const Real f = (q_point[gp](0) +
209  // q_point[gp](1) +
210  // q_point[gp](2));
211 
212  // add jac*weight*f*phi to the RHS in position ig
213 
214  Fu(i) += JxW[gp]*f*phi[i][gp];
215  Fv(i) += JxW[gp]*f*phi[i][gp];
216 
217  for (std::size_t j=0; j<phi.size(); ++j)
218  {
219 
220  Kuu(i,j) += JxW[gp]*((phi[i][gp])*(phi[j][gp]));
221 
222  Kvv(i,j) += JxW[gp]*((phi[i][gp])*(phi[j][gp]) +
223  1.*((dphi[i][gp])*(dphi[j][gp])));
224  };
225  };
226  vol += JxW[gp];
227  };
228 
229 
230  // You can't compute "area" (perimeter) if you are in 2D
231  if (dim == 3)
232  {
233  for (auto side : elem->side_index_range())
234  if (elem->neighbor_ptr(side) == libmesh_nullptr)
235  {
236  fe_face->reinit (elem, side);
237 
238  //fe_face->print_info();
239 
240  for (std::size_t gp=0; gp<JxW_face.size(); gp++)
241  area += JxW_face[gp];
242  }
243  }
244 
245  // Constrain the DOF indices.
246  dof_map.constrain_element_matrix_and_vector(Kuu, Fu, dof_indices_U);
247  dof_map.constrain_element_matrix_and_vector(Kvv, Fv, dof_indices_V);
248 
249 
250  system.rhs->add_vector(Fu, dof_indices_U);
251  system.rhs->add_vector(Fv, dof_indices_V);
252 
253  system.matrix->add_matrix(Kuu, dof_indices_U);
254  system.matrix->add_matrix(Kvv, dof_indices_V);
255  }
256 
257  libMesh::out << "Vol=" << vol << std::endl;
258 
259  if (dim == 3)
260  libMesh::out << "Area=" << area << std::endl;
261 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
unsigned int dim
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
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
void set_refinement_flag(const RefinementState rflag)
Sets the value of the refinement flag for the element.
Definition: elem.h:2513
NumericVector< Number > * rhs
The system matrix.
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 writing meshes in the GMV format.
Definition: gmv_io.h:46
This is the MeshBase class.
Definition: mesh_base.h:68
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
void assemble(EquationSystems &es, const std::string &system_name)
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.
void print_dof_constraints(std::ostream &os=libMesh::out, bool print_nonlocal=false) const
Prints (from processor 0) all DoF and Node constraints.
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
virtual void reinit()
Reinitialize all the systems.
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
CouplingMatrix * _dof_coupling
Degree of freedom coupling.
Definition: dof_map.h:1240
PetscErrorCode Vec Mat libmesh_dbg_var(j)
virtual void write(const std::string &) libmesh_override
This method implements writing a mesh to a specified file.
Definition: gmv_io.C:267
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 const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:490
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void resize(const unsigned int n)
Resizes the matrix and initializes all entries to be 0.
int main(int argc, char **argv)
Definition: amr.C:46
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
bool refine_and_coarsen_elements()
Refines and coarsens user-requested elements.
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.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
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.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
unsigned int n_points() const
Definition: quadrature.h:113
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.