GetFEM  5.4.2
helmholtz.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 /**
23  @file helmholtz.cc
24  @brief Helmholtz problem (Delta(u) + k^2 u = 0)
25 
26  Diffraction of a plane wave by a circular obstacle.
27 
28  This program is used to check that getfem++ is working. This is also
29  a good example of use of GetFEM.
30 */
31 
32 #include "getfem/getfem_assembling.h" /* import assembly methods (and norms comp.) */
33 #include "getfem/getfem_export.h" /* export functions (save solution in a file) */
36 #include "gmm/gmm.h"
37 using std::endl; using std::cout; using std::cerr;
38 using std::ends; using std::cin;
39 
40 /* some GetFEM types that we will be using */
41 using bgeot::base_small_vector; /* special class for small (dim<16) vectors */
42 using bgeot::base_node; /* geometrical nodes(derived from base_small_vector)*/
43 using bgeot::scalar_type; /* = double */
44 using bgeot::complex_type; /* = std::complex<double> */
45 using bgeot::size_type; /* = unsigned long */
46 using bgeot::short_type;
47 using bgeot::base_matrix; /* small dense matrix. */
48 
49 /* definition of some matrix/vector types. These ones are built
50  using the predefined types in Gmm++ */
52 typedef getfem::modeling_standard_complex_sparse_matrix sparse_matrix;
53 typedef getfem::modeling_standard_complex_plain_vector plain_vector;
54 
55 /*
56  structure for the Helmholtz problem
57 */
58 struct Helmholtz_problem {
59 
60  enum { DIRICHLET_BOUNDARY_NUM = 0, ROBIN_BOUNDARY_NUM = 1};
61  getfem::mesh mesh; /* the mesh */
62  getfem::mesh_im mim; /* the integration methods */
63  getfem::mesh_fem mf_u; /* main mesh_fem, for the elastostatic solution */
64  getfem::mesh_fem mf_rhs; /* mesh_fem for the right hand side (f(x),..) */
65  complex_type wave_number;
66 
67  scalar_type residual; /* max residual for the iterative solvers */
68  int with_mult;
69 
70  std::string datafilename;
71  bgeot::md_param PARAM;
72 
73  bool solve(plain_vector &U);
74  void init(void);
75  void compute_error(plain_vector &U);
76  Helmholtz_problem(void) : mim(mesh), mf_u(mesh), mf_rhs(mesh) {}
77 };
78 
79 complex_type __wave_number;
80 
81 complex_type incoming_field(const base_node &P) {
82  return complex_type(cos(__wave_number.real()*P[1]+.2),
83  sin(__wave_number.real()*P[1]+.2));
84  /*scalar_type s = 0;
85  for (size_type i=1; i < P.size(); ++i) s += P[i]*(1.-P[i]);
86  s = rand()*3. / RAND_MAX;
87  return s;
88  */
89 }
90 
91 /* Read parameters from the .param file, build the mesh, set finite element
92  * and integration methods and selects the boundaries.
93  */
94 void Helmholtz_problem::init(void) {
95  std::string FEM_TYPE = PARAM.string_value("FEM_TYPE","FEM name");
96  std::string INTEGRATION = PARAM.string_value("INTEGRATION",
97  "Name of integration method");
98  cout << "FEM_TYPE=" << FEM_TYPE << "\n";
99  cout << "INTEGRATION=" << INTEGRATION << "\n";
100 
101  /* First step : build the mesh */
102  size_type Nt = PARAM.int_value("NTHETA", "Nomber of space steps "),
103  Nr=PARAM.int_value("NR", "Nomber of space steps ");
104  size_type gt_order = PARAM.int_value("GTDEGREE",
105  "polynomial degree of geometric transformation");
106  scalar_type dtheta=2*M_PI*1./scalar_type(Nt);
107  scalar_type R0 = PARAM.real_value("R0","R0");
108  scalar_type R1 = PARAM.real_value("R1","R1");
109  scalar_type dR = (R1-R0)/scalar_type(Nr-1);
111  pgt = bgeot::parallelepiped_geotrans(2, short_type(gt_order));
112  for (size_type i=0; i < Nt; ++i) {
113  for (size_type j=0; j < Nr-1; ++j) {
114  std::vector<size_type> ipts; ipts.reserve(gmm::sqr(gt_order+1));
115  for (size_type ii=0; ii <= gt_order; ++ii) {
116  for (size_type jj=0; jj <= gt_order; ++jj) {
117  scalar_type r = R0 + scalar_type(j)*dR
118  + scalar_type(jj)*(dR/scalar_type(gt_order));
119  scalar_type t = scalar_type(i)*dtheta
120  + scalar_type(ii)*dtheta/scalar_type(gt_order);
121  ipts.push_back(mesh.add_point(base_node(r*cos(t),r*sin(t))));
122  }
123  }
124  mesh.add_convex(pgt, ipts.begin());
125  }
126  }
127 
128  datafilename = PARAM.string_value("ROOTFILENAME","Base name of data files.");
129  residual = PARAM.real_value("RESIDUAL"); if (residual == 0.) residual = 1e-10;
130 
131  __wave_number = wave_number = complex_type
132  (PARAM.real_value("WAVENUM_R", "Real part of the wave number"),
133  PARAM.real_value("WAVENUM_I", "Imaginary part of the wave number"));
134 
135  with_mult = int(PARAM.int_value("DIRICHLET_VERSION",
136  "Dirichlet condition version"));
137 
138  /* set the finite element on the mf_u */
139  getfem::pfem pf_u =
140  getfem::fem_descriptor(FEM_TYPE);
141  getfem::pintegration_method ppi =
142  getfem::int_method_descriptor(INTEGRATION);
143 
144  mim.set_integration_method(ppi);
145  mf_u.set_finite_element(pf_u);
146 
147  /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is
148  not used in the .param file */
149  std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE");
150  if (data_fem_name.size() == 0) {
151  GMM_ASSERT1(pf_u->is_lagrange(), "You are using a non-lagrange FEM. "
152  << "In that case you need to set "
153  << "DATA_FEM_TYPE in the .param file");
154  mf_rhs.set_finite_element(pf_u);
155  } else {
156  mf_rhs.set_finite_element(getfem::fem_descriptor(data_fem_name));
157  }
158 
159 
160  /* select boundaries */
161  cout << "Selecting Robin and Dirichlet boundaries\n";
162  getfem::mesh_region border_faces;
163  getfem::outer_faces_of_mesh(mesh, border_faces);
164  for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) {
165  assert(i.is_face());
166  if (gmm::vect_norm2(mesh.points_of_face_of_convex(i.cv(),
167  i.f())[0]) > 5.) {
168  mesh.region(ROBIN_BOUNDARY_NUM).add(i.cv(),i.f());
169  } else mesh.region(DIRICHLET_BOUNDARY_NUM).add(i.cv(),i.f());
170  }
171 }
172 
173 /**************************************************************************/
174 /* Model. */
175 /**************************************************************************/
176 
177 
178 bool Helmholtz_problem::solve(plain_vector &U) {
179 
180  // Complex model.
181  getfem::model Helmholtz_model(true);
182 
183  // Main unknown of the problem
184  Helmholtz_model.add_fem_variable("u", mf_u);
185 
186  // Helmholtz term on u.
187  plain_vector K(1); K[0] = wave_number;
188  Helmholtz_model.add_initialized_fixed_size_data("k", K);
189  add_Helmholtz_brick(Helmholtz_model, mim, "u", "k");
190 
191  // Fourier-Robin condition.
192  plain_vector Q(1); Q[0] = wave_number * complex_type(0,1.);
193  Helmholtz_model.add_initialized_fixed_size_data("Q", Q);
194  add_Fourier_Robin_brick(Helmholtz_model, mim, "u", "Q", ROBIN_BOUNDARY_NUM);
195 
196  // Dirichlet condition
197  plain_vector F(mf_rhs.nb_dof());
198  getfem::interpolation_function(mf_rhs, F, incoming_field);
199  Helmholtz_model.add_initialized_fem_data("DirichletData", mf_rhs, F);
201  (Helmholtz_model, mim, "u", mf_u,
202  DIRICHLET_BOUNDARY_NUM, "DirichletData");
203 
204  // Helmholtz_model.listvar(cout);
205 
206  gmm::iteration iter(residual, 1, 40000);
207  getfem::standard_solve(Helmholtz_model, iter);
208 
209  gmm::resize(U, mf_u.nb_dof());
210  gmm::copy(Helmholtz_model.complex_variable("u"), U);
211 
212  // cout << "U = " << U << endl;
213 
214  return (iter.converged());
215 }
216 
217 /**************************************************************************/
218 /* main program. */
219 /**************************************************************************/
220 
221 int main(int argc, char *argv[]) {
222 
223  GETFEM_MPI_INIT(argc, argv);
224  GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
225  FE_ENABLE_EXCEPT; // Enable floating point exception for Nan.
226 
227  Helmholtz_problem p;
228  p.PARAM.read_command_line(argc, argv);
229  p.init();
230  plain_vector U(p.mf_u.nb_dof());
231  if (!p.solve(U)) GMM_ASSERT1(false, "Solve has failed");
232 
233  if (p.PARAM.int_value("VTK_EXPORT")) {
234  cout << "export to " << p.datafilename + ".vtk" << "..\n";
235  getfem::vtk_export vtk_exp(p.datafilename + ".vtk",
236  p.PARAM.int_value("VTK_EXPORT")==1, true);
237  cout << "export to " << p.datafilename + ".vtu" << "..\n";
238  getfem::vtu_export vtu_exp(p.datafilename + ".vtu");
239  getfem::stored_mesh_slice sl(p.mesh, p.mesh.nb_convex() < 2000 ? 8 : 6);
240  vtk_exp.exporting(sl);
241  vtk_exp.write_point_data(p.mf_u, gmm::real_part(U), "helmholtz_rfield");
242  vtk_exp.write_point_data(p.mf_u, gmm::imag_part(U), "helmholtz_ifield");
243  vtu_exp.exporting(sl);
244  vtu_exp.write_point_data(p.mf_u, gmm::real_part(U), "helmholtz_rfield");
245  vtu_exp.write_point_data(p.mf_u, gmm::imag_part(U), "helmholtz_ifield");
246  cout << "export done, you can view the data file with (for example)\n"
247  "mayavi2 -d helmholtz.vtk -f WarpScalar -m Surface -m Outline"
248  "\n";
249  }
250 
251  GETFEM_MPI_FINALIZE;
252 
253  return 0;
254 }
getfem_export.h
Export solutions to various formats.
getfem_model_solvers.h
Standard solvers for model bricks.
getfem::stored_mesh_slice
The output of a getfem::mesh_slicer which has been recorded.
Definition: getfem_mesh_slice.h:47
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
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.
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::add_Dirichlet_condition_with_multipliers
size_type APIDECL add_Dirichlet_condition_with_multipliers(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region, const std::string &dataname=std::string())
Add a Dirichlet condition on the variable varname and the mesh region region.
Definition: getfem_models.cc:4671
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::add_Helmholtz_brick
size_type APIDECL add_Helmholtz_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region=size_type(-1))
Add a Helmoltz brick to the model.
Definition: getfem_models.cc:5516
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
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
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::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
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
getfem::standard_solve
void standard_solve(model &md, gmm::iteration &iter, rmodel_plsolver_type lsolver, abstract_newton_line_search &ls)
A default solver for the model brick system.
Definition: getfem_model_solvers.cc:410
getfem::add_Fourier_Robin_brick
size_type APIDECL add_Fourier_Robin_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region)
Add a Fourier-Robin brick to the model.
Definition: getfem_models.cc:5637
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....
getfem::vtk_export
VTK/VTU export.
Definition: getfem_export.h:68