1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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
main(int argc,char ** argv)46 int main (int argc, char ** argv)
47 {
48 LibMeshInit init(argc, argv);
49
50 libmesh_error_msg_if(argc < 4, "Usage: ./prog -d DIM filename");
51
52 // Variables to get us started
53 const unsigned char dim = cast_int<unsigned char>(atoi(argv[2]));
54
55 std::string meshname (argv[3]);
56
57 // declare a mesh...
58 Mesh mesh(init.comm(), dim);
59
60 // Read a mesh
61 mesh.read(meshname);
62
63 GMVIO(mesh).write ("out_0.gmv");
64
65 mesh.elem_ref(0).set_refinement_flag (Elem::REFINE);
66
67 MeshRefinement mesh_refinement (mesh);
68
69 mesh_refinement.refine_and_coarsen_elements ();
70 mesh_refinement.uniformly_refine (2);
71
72 mesh.print_info();
73
74
75 // Set up the equation system(s)
76 EquationSystems es (mesh);
77
78 LinearImplicitSystem & primary =
79 es.add_system<LinearImplicitSystem>("primary");
80
81 primary.add_variable ("U", FIRST);
82 primary.add_variable ("V", FIRST);
83
84 primary.get_dof_map()._dof_coupling->resize(2);
85 (*primary.get_dof_map()._dof_coupling)(0,0) = 1;
86 (*primary.get_dof_map()._dof_coupling)(1,1) = 1;
87
88 primary.attach_assemble_function(assemble);
89
90 es.init ();
91
92 es.print_info ();
93 primary.get_dof_map().print_dof_constraints ();
94
95 // call the solver.
96 primary.solve ();
97
98 GMVIO(mesh).write_equation_systems ("out_1.gmv",
99 es);
100
101
102
103 // Refine uniformly
104 mesh_refinement.uniformly_refine (1);
105 es.reinit ();
106
107 // Write out the projected solution
108 GMVIO(mesh).write_equation_systems ("out_2.gmv",
109 es);
110
111 // Solve again. Output the refined solution
112 primary.solve ();
113 GMVIO(mesh).write_equation_systems ("out_3.gmv",
114 es);
115
116 return 0;
117 }
118 #else
main(int,char **)119 int main (int, char **)
120 {
121 return 1;
122 }
123 #endif // ENABLE_AMR
124
125
126
127
128
129
assemble(EquationSystems & es,const std::string & libmesh_dbg_var (system_name))130 void assemble(EquationSystems & es,
131 const std::string & libmesh_dbg_var(system_name))
132 {
133 libmesh_assert_equal_to (system_name, "primary");
134
135 const MeshBase & mesh = es.get_mesh();
136 const unsigned int dim = mesh.mesh_dimension();
137
138 // Also use a 3x3x3 quadrature rule (3D). Then tell the FE
139 // about the geometry of the problem and the quadrature rule
140 FEType fe_type (FIRST);
141
142 std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
143 QGauss qrule(dim, FIFTH);
144
145 fe->attach_quadrature_rule (&qrule);
146
147 std::unique_ptr<FEBase> fe_face(FEBase::build(dim, fe_type));
148 QGauss qface(dim-1, FIFTH);
149
150 fe_face->attach_quadrature_rule(&qface);
151
152 LinearImplicitSystem & system =
153 es.get_system<LinearImplicitSystem>("primary");
154
155
156 // These are references to cell-specific data
157 const std::vector<Real> & JxW_face = fe_face->get_JxW();
158 const std::vector<Real> & JxW = fe->get_JxW();
159 const std::vector<Point> & q_point = fe->get_xyz();
160 const std::vector<std::vector<Real>> & phi = fe->get_phi();
161 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
162
163 std::vector<dof_id_type> dof_indices_U;
164 std::vector<dof_id_type> dof_indices_V;
165 const DofMap & dof_map = system.get_dof_map();
166
167 DenseMatrix<Number> Kuu;
168 DenseMatrix<Number> Kvv;
169 DenseVector<Number> Fu;
170 DenseVector<Number> Fv;
171
172 Real vol=0., area=0.;
173
174 for (const auto & elem : mesh.active_local_element_ptr_range())
175 {
176 // recompute the element-specific data for the current element
177 fe->reinit (elem);
178
179
180 //fe->print_info();
181
182 dof_map.dof_indices(elem, dof_indices_U, 0);
183 dof_map.dof_indices(elem, dof_indices_V, 1);
184
185 const unsigned int n_phi = cast_int<unsigned int>(phi.size());
186
187 // zero the element matrix and vector
188 Kuu.resize (n_phi, n_phi);
189
190 Kvv.resize (n_phi, n_phi);
191
192 Fu.resize (n_phi);
193 Fv.resize (n_phi);
194
195 // standard stuff... like in code 1.
196 for (unsigned int gp=0; gp<qrule.n_points(); gp++)
197 {
198 for (unsigned int i=0; i<n_phi; ++i)
199 {
200 // this is tricky. ig is the _global_ dof index corresponding
201 // to the _global_ vertex number elem->node_id(i). Note that
202 // in general these numbers will not be the same (except for
203 // the case of one unknown per node on one subdomain) so
204 // we need to go through the dof_map
205
206 const Real f = q_point[gp]*q_point[gp];
207 // const Real f = (q_point[gp](0) +
208 // q_point[gp](1) +
209 // q_point[gp](2));
210
211 // add jac*weight*f*phi to the RHS in position ig
212
213 Fu(i) += JxW[gp]*f*phi[i][gp];
214 Fv(i) += JxW[gp]*f*phi[i][gp];
215
216 for (unsigned int j=0; j != n_phi; ++j)
217 {
218
219 Kuu(i,j) += JxW[gp]*((phi[i][gp])*(phi[j][gp]));
220
221 Kvv(i,j) += JxW[gp]*((phi[i][gp])*(phi[j][gp]) +
222 1.*((dphi[i][gp])*(dphi[j][gp])));
223 };
224 };
225 vol += JxW[gp];
226 };
227
228
229 // You can't compute "area" (perimeter) if you are in 2D
230 if (dim == 3)
231 {
232 for (auto side : elem->side_index_range())
233 if (elem->neighbor_ptr(side) == nullptr)
234 {
235 fe_face->reinit (elem, side);
236
237 for (const auto & val : JxW_face)
238 area += val;
239 }
240 }
241
242 // Constrain the DOF indices.
243 dof_map.constrain_element_matrix_and_vector(Kuu, Fu, dof_indices_U);
244 dof_map.constrain_element_matrix_and_vector(Kvv, Fv, dof_indices_V);
245
246
247 system.rhs->add_vector(Fu, dof_indices_U);
248 system.rhs->add_vector(Fv, dof_indices_V);
249
250 system.matrix->add_matrix(Kuu, dof_indices_U);
251 system.matrix->add_matrix(Kvv, dof_indices_V);
252 }
253
254 libMesh::out << "Vol=" << vol << std::endl;
255
256 if (dim == 3)
257 libMesh::out << "Area=" << area << std::endl;
258 }
259