GetFEM  5.4.2
laplacian.cc
Go to the documentation of this file.
1 /*===========================================================================
2 
3  Copyright (C) 2002-2020 Yves Renard, Julien Pommier.
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 /**@file laplacian.cc
23  @brief Laplacian (Poisson) problem.
24 
25  The laplace equation is solved on a regular mesh of the unit
26  square, and is compared to the analytical solution.
27 
28  This program is used to check that getfem++ is working. This is
29  also a good example of use of GetFEM. This program does not use the
30  model bricks intentionally in order to serve as an example of solving
31  a pde directly with the assembly procedures.
32 */
33 
34 #include "getfem/getfem_assembling.h"
35 #include "getfem/getfem_export.h"
36 #include "getfem/getfem_regular_meshes.h"
37 #include "getfem/getfem_derivatives.h"
38 #include "getfem/getfem_superlu.h"
39 #include "gmm/gmm.h"
40 using std::endl; using std::cout; using std::cerr;
41 using std::ends; using std::cin;
42 
43 /* some GetFEM types that we will be using */
44 using bgeot::base_small_vector; /* special class for small (dim<16) vectors */
45 using bgeot::base_node; /* geometrical nodes (derived from base_small_vector)*/
46 using bgeot::scalar_type; /* = double */
47 using bgeot::size_type; /* = unsigned long */
48 
49 /* definition of some matrix/vector types. These ones are built
50  * using the predefined types in Gmm++
51  */
52 typedef gmm::rsvector<scalar_type> sparse_vector_type;
53 typedef gmm::row_matrix<sparse_vector_type> sparse_matrix_type;
54 typedef gmm::col_matrix<sparse_vector_type> col_sparse_matrix_type;
55 typedef std::vector<scalar_type> plain_vector;
56 
57 /* Definitions for the exact solution of the Laplacian problem,
58  * i.e. Delta(sol_u) + sol_f = 0
59  */
60 
61 base_small_vector sol_K; /* a coefficient */
62 /* exact solution */
63 scalar_type sol_u(const base_node &x) { return sin(gmm::vect_sp(sol_K, x)); }
64 /* righ hand side */
65 scalar_type sol_f(const base_node &x)
66 { return gmm::vect_sp(sol_K, sol_K) * sin(gmm::vect_sp(sol_K, x)); }
67 /* gradient of the exact solution */
68 base_small_vector sol_grad(const base_node &x)
69 { return sol_K * cos(gmm::vect_sp(sol_K, x)); }
70 
71 /*
72  structure for the Laplacian problem
73 */
74 struct laplacian_problem {
75 
76  enum { DIRICHLET_BOUNDARY_NUM = 0, NEUMANN_BOUNDARY_NUM = 1};
77  getfem::mesh mesh; /* the mesh */
78  getfem::mesh_im mim; /* the integration methods. */
79  getfem::mesh_fem mf_u; /* the main mesh_fem, for the Laplacian solution */
80  getfem::mesh_fem mf_rhs; /* the mesh_fem for the right hand side(f(x),..) */
81  getfem::mesh_fem mf_coef; /* the mesh_fem to represent pde coefficients */
82 
83  scalar_type residual; /* max residual for the iterative solvers */
84  size_type N;
85  bool gen_dirichlet;
86 
87  sparse_matrix_type SM; /* stiffness matrix. */
88  std::vector<scalar_type> U, B; /* main unknown, and right hand side */
89 
90  std::vector<scalar_type> Ud; /* reduced sol. for gen. Dirichlet condition. */
91  col_sparse_matrix_type NS; /* Dirichlet NullSpace
92  * (used if gen_dirichlet is true)
93  */
94  std::string datafilename;
95  bgeot::md_param PARAM;
96 
97  void assembly(void);
98  bool solve(void);
99  void init(void);
100  void compute_error();
101  laplacian_problem(void) : mim(mesh), mf_u(mesh), mf_rhs(mesh),
102  mf_coef(mesh) {}
103 };
104 
105 /* Read parameters from the .param file, build the mesh, set finite element
106  * and integration methods and selects the boundaries.
107  */
108 void laplacian_problem::init(void) {
109 
110  std::string MESH_TYPE = PARAM.string_value("MESH_TYPE","Mesh type ");
111  std::string FEM_TYPE = PARAM.string_value("FEM_TYPE","FEM name");
112  std::string INTEGRATION = PARAM.string_value("INTEGRATION",
113  "Name of integration method");
114 
115  cout << "MESH_TYPE=" << MESH_TYPE << "\n";
116  cout << "FEM_TYPE=" << FEM_TYPE << "\n";
117  cout << "INTEGRATION=" << INTEGRATION << "\n";
118 
119  /* First step : build the mesh */
120  bgeot::pgeometric_trans pgt =
121  bgeot::geometric_trans_descriptor(MESH_TYPE);
122  N = pgt->dim();
123  std::vector<size_type> nsubdiv(N);
124  std::fill(nsubdiv.begin(),nsubdiv.end(),
125  PARAM.int_value("NX", "Nomber of space steps "));
126  getfem::regular_unit_mesh(mesh, nsubdiv, pgt,
127  PARAM.int_value("MESH_NOISED") != 0);
128 
129  bgeot::base_matrix M(N,N);
130  for (size_type i=0; i < N; ++i) {
131  static const char *t[] = {"LX","LY","LZ"};
132  M(i,i) = (i<3) ? PARAM.real_value(t[i],t[i]) : 1.0;
133  }
134  if (N>1) { M(0,1) = PARAM.real_value("INCLINE") * PARAM.real_value("LY"); }
135 
136  /* scale the unit mesh to [LX,LY,..] and incline it */
137  mesh.transformation(M);
138 
139  datafilename = PARAM.string_value("ROOTFILENAME","Base name of data files.");
140  scalar_type FT = PARAM.real_value("FT", "parameter for exact solution");
141  residual = PARAM.real_value("RESIDUAL");
142  if (residual == 0.) residual = 1e-10;
143  sol_K.resize(N);
144  for (size_type j = 0; j < N; j++)
145  sol_K[j] = ((j & 1) == 0) ? FT : -FT;
146 
147  /* set the finite element on the mf_u */
148  getfem::pfem pf_u = getfem::fem_descriptor(FEM_TYPE);
149  getfem::pintegration_method ppi = getfem::int_method_descriptor(INTEGRATION);
150 
151  mim.set_integration_method(mesh.convex_index(), ppi);
152  mf_u.set_finite_element(mesh.convex_index(), pf_u);
153 
154  /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is
155  not used in the .param file */
156  std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE");
157  if (data_fem_name.size() == 0) {
158  GMM_ASSERT1(pf_u->is_lagrange(), "You are using a non-lagrange FEM. "
159  << "In that case you need to set "
160  << "DATA_FEM_TYPE in the .param file");
161  mf_rhs.set_finite_element(mesh.convex_index(), pf_u);
162  } else {
163  mf_rhs.set_finite_element(mesh.convex_index(),
164  getfem::fem_descriptor(data_fem_name));
165  }
166 
167  /* set the finite element on mf_coef. Here we use a very simple element
168  * since the only function that need to be interpolated on the mesh_fem
169  * is f(x)=1 ... */
170  mf_coef.set_finite_element(mesh.convex_index(),
171  getfem::classical_fem(pgt,0));
172 
173  /* set boundary conditions
174  * (Neuman on the upper face, Dirichlet elsewhere) */
175  gen_dirichlet = PARAM.int_value("GENERIC_DIRICHLET");
176 
177  if (!pf_u->is_lagrange() && !gen_dirichlet)
178  GMM_WARNING2("With non lagrange fem prefer the generic "
179  "Dirichlet condition option");
180 
181  cout << "Selecting Neumann and Dirichlet boundaries\n";
182  getfem::mesh_region border_faces;
183  getfem::outer_faces_of_mesh(mesh, border_faces);
184  for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) {
185  assert(i.is_face());
186  base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f());
187  un /= gmm::vect_norm2(un);
188  if (gmm::abs(un[N-1] - 1.0) < 1.0E-7) { // new Neumann face
189  mesh.region(NEUMANN_BOUNDARY_NUM).add(i.cv(), i.f());
190  } else {
191  mesh.region(DIRICHLET_BOUNDARY_NUM).add(i.cv(), i.f());
192  }
193  }
194 }
195 
196 void laplacian_problem::assembly(void) {
197  size_type nb_dof = mf_u.nb_dof();
198  size_type nb_dof_rhs = mf_rhs.nb_dof();
199 
200  gmm::resize(B, nb_dof); gmm::clear(B);
201  gmm::resize(U, nb_dof); gmm::clear(U);
202  gmm::resize(SM, nb_dof, nb_dof); gmm::clear(SM);
203 
204  cout << "Number of dof : " << nb_dof << endl;
205  cout << "Assembling stiffness matrix" << endl;
206  getfem::asm_stiffness_matrix_for_laplacian(SM, mim, mf_u, mf_coef,
207  std::vector<scalar_type>(mf_coef.nb_dof(), 1.0));
208 
209  cout << "Assembling source term" << endl;
210  std::vector<scalar_type> F(nb_dof_rhs);
211  getfem::interpolation_function(mf_rhs, F, sol_f);
212  getfem::asm_source_term(B, mim, mf_u, mf_rhs, F);
213 
214  cout << "Assembling Neumann condition" << endl;
215  gmm::resize(F, nb_dof_rhs*N);
216  getfem::interpolation_function(mf_rhs, F, sol_grad);
217  getfem::asm_normal_source_term(B, mim, mf_u, mf_rhs, F,
218  NEUMANN_BOUNDARY_NUM);
219 
220  cout << "take Dirichlet condition into account" << endl;
221  if (!gen_dirichlet) {
222  std::vector<scalar_type> D(nb_dof);
223  getfem::interpolation_function(mf_u, D, sol_u);
224  getfem::assembling_Dirichlet_condition(SM, B, mf_u,
225  DIRICHLET_BOUNDARY_NUM, D);
226  } else {
227  gmm::resize(F, nb_dof_rhs);
228  getfem::interpolation_function(mf_rhs, F, sol_u);
229 
230  gmm::resize(Ud, nb_dof);
231  gmm::resize(NS, nb_dof, nb_dof);
232  col_sparse_matrix_type H(nb_dof_rhs, nb_dof);
233  std::vector<scalar_type> R(nb_dof_rhs);
234  std::vector<scalar_type> RHaux(nb_dof);
235 
236  /* build H and R such that U mush satisfy H*U = R */
237  getfem::asm_dirichlet_constraints(H, R, mim, mf_u, mf_rhs,
238  mf_rhs, F, DIRICHLET_BOUNDARY_NUM);
239 
240  gmm::clean(H, 1e-12);
241 // cout << "H = " << H << endl;
242 // cout << "R = " << R << endl;
243  int nbcols = int(getfem::Dirichlet_nullspace(H, NS, R, Ud));
244  // cout << "Number of irreductible unknowns : " << nbcols << endl;
245  gmm::resize(NS, gmm::mat_ncols(H),nbcols);
246 
247  gmm::mult(SM, Ud, gmm::scaled(B, -1.0), RHaux);
248  gmm::resize(B, nbcols);
249  gmm::resize(U, nbcols);
250  gmm::mult(gmm::transposed(NS), gmm::scaled(RHaux, -1.0), B);
251  sparse_matrix_type SMaux(nbcols, nb_dof);
252  gmm::mult(gmm::transposed(NS), SM, SMaux);
253  gmm::resize(SM, nbcols, nbcols);
254  /* NSaux = NS, but is stored by rows instead of by columns */
255  sparse_matrix_type NSaux(nb_dof, nbcols); gmm::copy(NS, NSaux);
256  gmm::mult(SMaux, NSaux, SM);
257  }
258 }
259 
260 
261 bool laplacian_problem::solve(void) {
262 
263  // see_schmidt(SM, U, B);
264 
265  cout << "Compute preconditionner\n";
266  gmm::iteration iter(residual, 1, 40000);
267  double time = gmm::uclock_sec();
268  if (1) {
269  // gmm::identity_matrix P;
270  // gmm::diagonal_precond<sparse_matrix_type> P(SM);
271  // gmm::mr_approx_inverse_precond<sparse_matrix_type> P(SM, 10, 10E-17);
272  // gmm::ildlt_precond<sparse_matrix_type> P(SM);
273  // gmm::ildltt_precond<sparse_matrix_type> P(SM, 20, 1E-6);
274  gmm::ilut_precond<sparse_matrix_type> P(SM, 20, 1E-6);
275  // gmm::ilutp_precond<sparse_matrix_type> P(SM, 20, 1E-6);
276  // gmm::ilu_precond<sparse_matrix_type> P(SM);
277  cout << "Time to compute preconditionner : "
278  << gmm::uclock_sec() - time << " seconds\n";
279 
280 
281  //gmm::HarwellBoeing_IO::write("SM", SM);
282 
283  // gmm::cg(SM, U, B, P, iter);
284  gmm::gmres(SM, U, B, P, 50, iter);
285  } else {
286  double rcond;
287  gmm::SuperLU_solve(SM, U, B, rcond);
288  cout << "cond = " << 1/rcond << "\n";
289  }
290 
291  cout << "Total time to solve : "
292  << gmm::uclock_sec() - time << " seconds\n";
293 
294  if (gen_dirichlet) {
295  std::vector<scalar_type> Uaux(mf_u.nb_dof());
296  gmm::mult(NS, U, Ud, Uaux);
297  gmm::resize(U, mf_u.nb_dof());
298  gmm::copy(Uaux, U);
299  }
300 
301  return (iter.converged());
302 }
303 
304 /* compute the error with respect to the exact solution */
305 void laplacian_problem::compute_error() {
306  std::vector<scalar_type> V(mf_rhs.nb_basic_dof());
307  getfem::interpolation(mf_u, mf_rhs, U, V);
308  for (size_type i = 0; i < mf_rhs.nb_basic_dof(); ++i)
309  V[i] -= sol_u(mf_rhs.point_of_basic_dof(i));
310  cout.precision(16);
311  cout << "L2 error = " << getfem::asm_L2_norm(mim, mf_rhs, V) << endl
312  << "H1 error = " << getfem::asm_H1_norm(mim, mf_rhs, V) << endl
313  << "Linfty error = " << gmm::vect_norminf(V) << endl;
314 }
315 
316 /**************************************************************************/
317 /* main program. */
318 /**************************************************************************/
319 
320 int main(int argc, char *argv[]) {
321 
322  GETFEM_MPI_INIT(argc, argv);
323  GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
324  FE_ENABLE_EXCEPT; // Enable floating point exception for Nan.
325 
326  try {
327  laplacian_problem p;
328  p.PARAM.read_command_line(argc, argv);
329  p.init();
330  p.mesh.write_to_file(p.datafilename + ".mesh");
331  p.assembly();
332  if (!p.solve()) GMM_ASSERT1(false, "Solve procedure has failed");
333  p.compute_error();
334  }
335  GMM_STANDARD_CATCH_ERROR;
336 
337  GETFEM_MPI_FINALIZE;
338 
339  return 0;
340 }
getfem_export.h
Export solutions to various formats.
getfem::asm_H1_norm
scalar_type asm_H1_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute the H1 norm of U.
Definition: getfem_assembling.h:302
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
getfem::asm_stiffness_matrix_for_laplacian
void asm_stiffness_matrix_for_laplacian(MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
assembly of , where is scalar.
Definition: getfem_assembling.h:1152
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::int_method_descriptor
pintegration_method int_method_descriptor(std::string name, bool throw_if_not_found=true)
Get an integration method from its name .
Definition: getfem_integration.cc:1130
getfem_regular_meshes.h
Build regular meshes.
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::Dirichlet_nullspace
size_type Dirichlet_nullspace(const MAT1 &H, MAT2 &NS, const VECT1 &R, VECT2 &U0)
Build an orthogonal basis of the kernel of H in NS, gives the solution of minimal norm of H*U = R in ...
Definition: getfem_assembling.h:1656
bgeot::geometric_trans_descriptor
pgeometric_trans geometric_trans_descriptor(std::string name)
Get the geometric transformation from its string name.
Definition: bgeot_geometric_trans.cc:1163
getfem::classical_fem
pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k, bool complete=false)
Give a pointer on the structures describing the classical polynomial fem of degree k on a given conve...
Definition: getfem_fem.cc:4141
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::interpolation
void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target, const VECTU &U, VECTV &V, int extrapolation=0, double EPS=1E-10, mesh_region rg_source=mesh_region::all_convexes(), mesh_region rg_target=mesh_region::all_convexes())
interpolation/extrapolation of (mf_source, U) on mf_target.
Definition: getfem_interpolation.h:693
getfem_superlu.h
SuperLU interface for getfem.
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
getfem::asm_normal_source_term
void asm_normal_source_term(VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg)
Normal source term (for boundary (Neumann) condition).
Definition: getfem_assembling.h:905
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
getfem::asm_source_term
void asm_source_term(const VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg=mesh_region::all_convexes())
source term (for both volumic sources and boundary (Neumann) sources).
Definition: getfem_assembling.h:877
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
getfem::asm_dirichlet_constraints
void asm_dirichlet_constraints(MAT &H, VECT1 &R, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_mult, const mesh_fem &mf_r, const VECT2 &r_data, const mesh_region &region, int version=ASMDIR_BUILDALL)
Assembly of Dirichlet constraints in a weak form.
Definition: getfem_assembling.h:1373
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
bgeot::small_vector
container for small vectors of POD (Plain Old Data) types.
Definition: bgeot_small_vector.h:205
getfem_derivatives.h
Compute the gradient of a field on a getfem::mesh_fem.
getfem::asm_L2_norm
scalar_type asm_L2_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute , U might be real or complex
Definition: getfem_assembling.h:56
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_def.h:488
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
gmm.h
Include common gmm files.
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::regular_unit_mesh
void regular_unit_mesh(mesh &m, std::vector< size_type > nsubdiv, bgeot::pgeometric_trans pgt, bool noised=false)
Build a regular mesh of the unit square/cube/, etc.
Definition: getfem_regular_meshes.cc:238
getfem::outer_faces_of_mesh
void APIDECL outer_faces_of_mesh(const mesh &m, const dal::bit_vector &cvlst, convex_face_ct &flist)
returns a list of "exterior" faces of a mesh (i.e.
Definition: getfem_mesh.cc:822
gmm::vect_norminf
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norminf(const V &v)
Infinity norm of a vector.
Definition: gmm_blas.h:693
getfem::interpolation_function
void interpolation_function(mesh_fem &mf_target, const VECT &VV, F &f, mesh_region rg=mesh_region::all_convexes())
interpolation of a function f on mf_target.
Definition: getfem_interpolation.h:186
gmm::ilut_precond
Incomplete LU with threshold and K fill-in Preconditioner.
Definition: gmm_precond_ilut.h:102
getfem::fem_descriptor
pfem fem_descriptor(const std::string &name)
get a fem descriptor from its string name.
Definition: getfem_fem.cc:4232
getfem_assembling.h
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....