Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
meshadaptation.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Cecile Daversin <cecile.daversin@lncmi.cnrs.fr>
6  Date: 2011-16-12
7 
8  Copyright (C) 2008-2010 Universite Joseph Fourier (Grenoble I)
9  Copyright (C) CNRS 2012
10 
11  This program is free software: you can redistribute it and/or modify
12  it under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  This program is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  GNU General Public License for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with this program. If not, see <http://www.gnu.org/licenses/>.
23 */
30 #ifndef __MESHADAPTATION_HPP
31 #define __MESHADAPTATION_HPP 1
32 
34 #include <feel/feelalg/backend.hpp>
35 
37 #include <feel/feeldiscr/functionspace.hpp>
38 
41 
43 #include <feel/feelvf/vf.hpp>
44 
45 #include <feel/feelcore/feel.hpp>
46 #include <feel/feeldiscr/mesh.hpp>
47 #include <feel/feeldiscr/projector.hpp>
48 
49 #include <boost/assert.hpp>
50 #include <boost/regex.hpp>
51 
52 #include <fstream>
53 
54 #if defined( FEELPP_HAS_GMSH_H )
55 #include <Gmsh.h>
56 #include <GModel.h>
57 #include <Context.h>
58 #include <MElement.h>
59 #include <OpenFile.h>
60 #include <StringUtils.h>
61 #include <Field.h>
62 #include <PView.h>
63 #include <PViewData.h>
64 #endif
65 
66 namespace Feel
67 {
68 
69  BOOST_PARAMETER_NAME(initMesh)
70  BOOST_PARAMETER_NAME(geofile)
71  BOOST_PARAMETER_NAME(adaptType)
72  BOOST_PARAMETER_NAME(var)
73  BOOST_PARAMETER_NAME(metric)
74  BOOST_PARAMETER_NAME(hMin)
75  BOOST_PARAMETER_NAME(hMax)
76  BOOST_PARAMETER_NAME(tol)
77 
78  template<int Dim,
79  int Order,
80  int OrderGeo>
81  class MeshAdaptation
82  {
83  public:
84  static const bool isP1 = (Order == 1);
85 
86  //typedef double value_type;
87  typedef Eigen::Matrix<double, Dim, 1> vectorN_type;
88  typedef Eigen::Matrix<double, Dim, Dim> matrixN_type;
89 
91  typedef Backend<double> backend_type;
93  typedef boost::shared_ptr<backend_type> backend_ptrtype;
94 
96  typedef typename backend_type::sparse_matrix_type sparse_matrix_type;
98  typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
100  typedef typename backend_type::vector_type vector_type;
102  typedef typename backend_type::vector_ptrtype vector_ptrtype;
103 
104  typedef Simplex<Dim,OrderGeo> convex_type;
106  typedef Mesh<convex_type> mesh_type;
108  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
109 
111  typedef Exporter<mesh_type> export_type;
113  typedef boost::shared_ptr<export_type> export_ptrtype;
114 
115  typedef bases<Lagrange<Order, Scalar> > basis_type;
116  typedef FunctionSpace<mesh_type, basis_type> space_type;
118  typedef boost::shared_ptr<space_type> space_ptrtype;
120  typedef typename space_type::element_type element_type;
121 
123  typedef bases<Lagrange<0,Scalar, Discontinuous > > p0_basis_type;
124  typedef FunctionSpace<mesh_type, p0_basis_type> p0_space_type;
125  typedef boost::shared_ptr<p0_space_type> p0_space_ptrtype;
126  typedef typename p0_space_type::element_type p0_element_type;
127 
129  typedef bases<Lagrange<1, Scalar> > p1_basis_type;
130  typedef FunctionSpace<mesh_type, p1_basis_type> p1_space_type;
131  typedef boost::shared_ptr<p1_space_type> p1_space_ptrtype;
132  typedef typename p1_space_type::element_type p1_element_type;
133 
135  typedef bases<Lagrange<1, Vectorial> > p1vec_basis_type;
136  typedef FunctionSpace<mesh_type, p1vec_basis_type> p1vec_space_type;
137  typedef boost::shared_ptr<p1vec_space_type> p1vec_space_ptrtype;
138  typedef typename p1vec_space_type::element_type p1vec_element_type;
139 
140  typedef typename mesh_type::point_type point_type;
141 
143  MeshAdaptation()
144  {
145  // build default value for _var parameter of interface
146  element_type initVar;
147  std::pair<element_type, std::string> initPairVar = std::make_pair( initVar, "defaultVar");
148  M_defaultVar.push_back(initPairVar);
149 
150  // build default value for _metric parameter of interface
151  p1_element_type initMetric1;
152  std::vector<p1_element_type> initMetric;
153  initMetric.push_back(initMetric1);
154  std::pair<std::vector<p1_element_type>, std::string> initPairMetric = std::make_pair(initMetric, "defaultMetric");
155  M_defaultMetric.push_back(initPairMetric);
156  }
157 
159  std::vector<vectorN_type> const measures()
160  {
161  return M_measures;
162  }
163 
165  std::vector<matrixN_type> const directions()
166  {
167  return M_directions;
168  }
169 
171  mesh_ptrtype adaptMeshImpl( const mesh_ptrtype& initMesh, std::string geofile, std::string adaptType,
172  std::list< std::pair<element_type, std::string> > var,
173  std::list< std::pair< std::vector<p1_element_type>, std::string> > metric,
174  double hMin, double hMax, double tol);
175 
177  std::string createPosfile(std::string name_var, const p1_element_type& bbNewMap, const mesh_ptrtype& mesh);
178  std::string createPosfileAnisotropic(std::string nameVar, const std::vector<p1_element_type>& bbNewMap, const mesh_ptrtype& mesh);
179 
181  std::string createAdaptedGeo(std::string geofile, std::string name, std::vector<std::string> posfiles, bool aniso);
182 
184  std::string buildAdaptedMesh(std::string geofile, std::string name, std::vector<std::string> posfiles , bool aniso);
185 
187  void computeMetric(const double tol, const double h_min, const double h_max, const matrixN_type & hessian_matrix,
188  matrixN_type & M, double & max_eigenvalue, int dofId);
189 
191  //1 : proj u Pk -> P1 and Hessian P1
192  //2 : Hessian P(k-2) and proj hessian P(k-2) -> P1 (possibility to choose L2 or H1 projection)
193  std::string adaptMeshHess1(element_type& U, const mesh_ptrtype& mesh, double hMin, double hMax,
194  std::string name, std::string geofile, double tol, bool aniso);
195 
196  std::string adaptMeshHess2(element_type& U, const mesh_ptrtype& mesh, double hMin, double hMax,
197  std::string name, std::string geofile, double tol, bool aniso);
198  std::string adaptMeshHess2(element_type& U, const mesh_ptrtype& mesh, double hMin, double hMax,
199  std::string name, std::string geofile, double tol, bool aniso,
200  mpl::bool_<true>);
201  std::string adaptMeshHess2(element_type& U, const mesh_ptrtype& mesh, double hMin, double hMax,
202  std::string name, std::string geofile, double tol, bool aniso,
203  mpl::bool_<false>);
204 
206  // _initMesh = initial mesh
207  // _geofile = complete path to access geofile
208  // _adaptType = isotropic | anisotropic
209  // _var = list of (element_type, string) = list of (variable, name of the variable)
210  // _metric = list of (vector<p1_element_type>, string) = list of (hsize(isotropic) | matrix(anisotropic), name for metric)
211  // _hMin, _hMax = limits for hsize
212  // _tol = geometric tolerance (for adaptation from hessian)
213 
214  BOOST_PARAMETER_MEMBER_FUNCTION(
215  (mesh_ptrtype),
216  adaptMesh,
217  tag,
218  (required
219  (initMesh,(mesh_ptrtype)) // initial mesh
220  (geofile, (std::string)) // geometry
221  (adaptType, (std::string))) //type of adaptation : isotropic | anisotropic
222  (optional
223  (var, (std::list< std::pair<element_type, std::string> >), M_defaultVar)
224  (metric, (std::list< std::pair< std::vector<p1_element_type>, std::string> >), M_defaultMetric)
225  (hMin, (double), 1.0e-3)
226  (hMax, (double), 100.0)
227  (tol, (double), 1.0))
228  )
229  {
230  mesh_ptrtype newMesh = this-> adaptMeshImpl( initMesh, geofile, adaptType, var, metric, hMin, hMax, tol);
231  return newMesh;
232  }
233 
234  private :
235  std::list< std::pair<element_type, std::string> > M_defaultVar;
236  std::list< std::pair<std::vector<p1_element_type>, std::string> > M_defaultMetric;
237 
238  std::vector<vectorN_type> M_measures;
239  std::vector<matrixN_type> M_directions;
240 
241  };
242 
243  template<int Dim,
244  int Order,
245  int OrderGeo>
246  const bool
247  MeshAdaptation<Dim, Order, OrderGeo>::isP1;
248 
249  template<int Dim,
250  int Order,
251  int OrderGeo>
252  typename MeshAdaptation<Dim, Order, OrderGeo>::mesh_ptrtype
253  MeshAdaptation<Dim,
254  Order,
255  OrderGeo>::adaptMeshImpl( const mesh_ptrtype& initMesh, std::string geofile, std::string adaptType,
256  std::list< std::pair<element_type, std::string> > var,
257  std::list< std::pair< std::vector<p1_element_type>, std::string> > metric,
258  double hMin, double hMax, double tol)
259  {
260  std::vector<std::string> posfiles;
261  std::string finalName = "";
262  if ( !(var.front().second == "defaultVar") ) // adptation from a variable
263  {
264  boost::for_each( var, [&initMesh, &hMin, &hMax, &geofile, &tol, &adaptType, &posfiles, &finalName, this]
265  (std::pair<element_type, std::string> _var)
266  {
267  std::string newPosfile = this->adaptMeshHess2( _var.first, initMesh, hMin, hMax, _var.second, geofile,
268  tol, adaptType=="anisotropic");
269  posfiles.push_back( newPosfile );
270  finalName += "-" + _var.second;
271  });
272  }
273 
274  else if ( !(metric.front().second == "defaultMetric") ) //adaptation from a metric
275  {
276  boost::for_each(metric, [&initMesh, &adaptType, &posfiles, &finalName, this]
277  (std::pair<std::vector<p1_element_type>, std::string> _metric)
278  {
279  std::string newPosfile;
280  if (adaptType == "isotropic")
281  {
282  BOOST_ASSERT(_metric.first.size()==1);
283  newPosfile = this->createPosfile(_metric.second, _metric.first[0], initMesh);
284  }
285  else if( adaptType == "anisotropic")
286  {
287  BOOST_ASSERT(_metric.first.size()==Dim*Dim);
288  newPosfile = this->createPosfileAnisotropic(_metric.second, _metric.first, initMesh);
289  }
290  posfiles.push_back( newPosfile );
291  finalName += "-" + _metric.second;
292  });
293  }
294 
295  std::string newMeshName = this->buildAdaptedMesh( geofile, finalName, posfiles, adaptType=="anisotropic");
296 
297  // Update the mesh
298  LOG(INFO) << "load the new mesh...\n";
299  mesh_ptrtype newMesh = mesh_type::New();
300  newMesh = loadGMSHMesh( _mesh = new mesh_type,
301  _filename = newMeshName,
302  _update=MESH_UPDATE_FACES | MESH_UPDATE_EDGES);
303  LOG(INFO) << "new mesh loaded \n";
304 
305  return newMesh;
306  }
307 
308  template<int Dim,
309  int Order,
310  int OrderGeo>
311  std::string
312  MeshAdaptation<Dim,
313  Order,
314  OrderGeo>::createPosfile(std::string nameVar, const p1_element_type& bbNewMap, const mesh_ptrtype& mesh)
315  {
316  using namespace Feel::vf;
317 
318  std::string posFormat = "pos";
319  std::string mshFormat = "msh";
320 
321  p1_space_ptrtype P1h = p1_space_type::New( mesh );
322 
323  std::ofstream newPosFile( (boost::format( "%1%.%2%" ) % nameVar %posFormat ).str() );
324 
325  newPosFile << "View \" background mesh \" { \n";
326  auto eltIt = mesh->beginElement();
327  auto eltEnd = mesh->endElement();
328  for ( ; eltIt != eltEnd; eltIt++)
329  {
330  std::vector<point_type> eltPoints( eltIt->nPoints() );
331  if (Dim == 2)
332  newPosFile << "ST(";
333  if (Dim == 3)
334  newPosFile << "SS(";
335  for (int i=0; i < eltIt->nPoints(); i++)
336  {
337  eltPoints[i] = eltIt->point(i);
338  for (int j=0; j< Dim; j++)
339  {
340  newPosFile << eltPoints[i](j);
341  if ( Dim == 2 || (Dim == 3 && j!=Dim-1) )
342  {
343  newPosFile << ", ";
344  }
345  }
346 
347  if (Dim == 2)
348  newPosFile << "0"; //2D case => z coordinate = 0
349  if (i!= eltIt->nPoints() - 1)
350  newPosFile << ", ";
351  }
352 
353  newPosFile << "){";
354 
355  for (size_t k=0; k<eltPoints.size(); k++)
356  {
357  auto dofIndex = boost::get<0>( P1h->dof()->localToGlobal( eltIt->id(), k) );
358  newPosFile << bbNewMap[ dofIndex ];
359  if ( k!= eltPoints.size() - 1)
360  newPosFile << ", ";
361  }
362 
363  newPosFile << "};\n";
364  eltPoints.clear();
365  }
366 
367  newPosFile << "};";
368  newPosFile.close();
369 
370  return (boost::format( "%1%.%2%" ) % nameVar %posFormat ).str();
371  }
372 
373  template<int Dim,
374  int Order,
375  int OrderGeo>
376  std::string
377  MeshAdaptation<Dim,
378  Order,
379  OrderGeo>::createPosfileAnisotropic(std::string nameVar, const std::vector<p1_element_type>& bbNewMap,
380  const mesh_ptrtype& mesh)
381  {
382  using namespace Feel::vf;
383 
384  std::string posFormat = "pos";
385  std::string mshFormat = "msh";
386 
387  p1_space_ptrtype P1h = p1_space_type::New( mesh );
388 
389  std::ofstream newPosFile( (boost::format( "%1%.%2%" ) % nameVar %posFormat ).str() );
390 
391  newPosFile << "View \" background mesh \" { \n";
392  auto eltIt = mesh->beginElement();
393  auto eltEnd = mesh->endElement();
394  for ( ; eltIt != eltEnd; eltIt++)
395  {
396  std::vector<point_type> eltPoints( eltIt->nPoints() );
397  newPosFile << "T";
398  if (Dim == 2)
399  newPosFile << "T";
400  else
401  newPosFile << "S";
402 
403  newPosFile << "(";
404  for (int i=0; i < eltIt->nPoints(); i++)
405  {
406  eltPoints[i] = eltIt->point(i);
407  for (int j=0; j< Dim; j++)
408  {
409  newPosFile << eltPoints[i](j);
410  if ( Dim == 2 || (Dim == 3 && j!=Dim-1) )
411  newPosFile << ", ";
412  }
413  if (Dim == 2)
414  newPosFile << "0"; //2D case => z coordinate = 0
415  if (i!= eltIt->nPoints() - 1)
416  newPosFile << ", ";
417  }
418 
419  newPosFile << "){";
420 
421  for (size_t k=0; k<eltPoints.size(); k++)
422  {
423  auto dofIndex = boost::get<0>( P1h->dof()->localToGlobal( eltIt->id(), k) );
424 
425  int num = 0;
426  newPosFile << (bbNewMap[num++])[ dofIndex ];
427  newPosFile << ", ";
428 
429  newPosFile << (bbNewMap[num++])[ dofIndex ] << ", ";
430 
431  if (Dim == 3)
432  newPosFile << (bbNewMap[num++])[ dofIndex ] << ", ";
433  else
434  newPosFile << "0, ";
435 
436  newPosFile << (bbNewMap[num++])[ dofIndex ] << ", ";
437  newPosFile << (bbNewMap[num++])[ dofIndex ] << ", ";
438 
439  if (Dim == 3)
440  {
441  newPosFile << (bbNewMap[num++])[ dofIndex ] << ", ";
442  newPosFile << (bbNewMap[num++])[ dofIndex ] << ", ";
443  newPosFile << (bbNewMap[num++])[ dofIndex ] << ", ";
444  newPosFile << (bbNewMap[num++])[ dofIndex ];
445  }
446  else
447  {
448  newPosFile << "0, ";
449  newPosFile << "0, ";
450  newPosFile << "0, ";
451  newPosFile << "1";
452  }
453 
454  if (k!= eltPoints.size() - 1)
455  newPosFile << ", ";
456  }
457  newPosFile << "};\n";
458  eltPoints.clear();
459  }
460 
461  newPosFile << "};";
462  newPosFile.close();
463 
464  return (boost::format( "%1%.%2%" ) % nameVar %posFormat ).str();
465  }
466 
467  template<int Dim,
468  int Order,
469  int OrderGeo>
470  std::string
471  MeshAdaptation<Dim,
472  Order,
473  OrderGeo>::createAdaptedGeo(std::string geofile, std::string name, std::vector<std::string> posfiles,
474  bool aniso)
475  {
477  std::ifstream stringToGeo;
478  std::string geofileString;
479  stringToGeo.open(geofile);
480  stringToGeo.seekg(0, std::ios::end);
481  geofileString.reserve(stringToGeo.tellg());
482  stringToGeo.seekg(0, std::ios::beg);
483  geofileString.assign((std::istreambuf_iterator<char>(stringToGeo)),
484  std::istreambuf_iterator<char>());
485  stringToGeo.close();
486 
487  int nbPosfiles = posfiles.size();
488 
490  std::ostringstream __exprAlgo, __exprFlist, __exprBgm;
491  __exprAlgo << "Field.([[:blank:]]*)([[:digit:]]*)([[:blank:]]*).([[:blank:]]*)=([[:blank:]]*)(Min|MinAniso);";
492  __exprFlist << "Field.([[:blank:]]*)([[:digit:]]*)([[:blank:]]*)..FieldsList([[:blank:]]*)=([[:blank:]]*).(([[:digit:]].?)*).;";
493  __exprBgm << "Background[[:blank:]]Field([[:blank:]]*)=([[:blank:]]*)([[:digit:]]*);";
494 
495  boost::regex exprAlgo(__exprAlgo.str().c_str());
496  boost::regex exprFlist(__exprFlist.str().c_str());
497  boost::regex exprBgm(__exprBgm.str().c_str());
498  std::list<boost::regex> endExpr = boost::assign::list_of(exprAlgo)(exprFlist)(exprBgm);
499 
500  bool isInserted = false;
501 
503  int i=1;
504  for (; i <= nbPosfiles ; i++)
505  {
506  std::ostringstream __expr1, __expr2;
507  __expr1 << "Field.([[:blank:]]*)"<< i <<"([[:blank:]]*).([[:blank:]]*)=([[:blank:]]*)PostView;" ;
508  __expr2 << "Field.([[:blank:]]*)" << i << "([[:blank:]]*)..IView([[:blank:]]*)=([[:blank:]]*)" << i-1 << ";" ;
509 
510  boost::regex expression1( __expr1.str().c_str() );
511  boost::regex expression2( __expr2.str().c_str() );
512 
513  boost::match_results<std::string::iterator> posFields;
514  bool fieldIsFoundL1 = boost::regex_search(geofileString.begin(), geofileString.end(), posFields, expression1);
515  bool fieldIsFoundL2 = boost::regex_search(geofileString.begin(), geofileString.end(), posFields, expression2);
516 
518  if (!(fieldIsFoundL1 && fieldIsFoundL2))
519  {
520  std::ostringstream __newExpr;
521  __newExpr << "Field[" << i << "] = PostView; \n";
522  __newExpr << "Field[" << i << "].IView = " << i-1 << ";\n";
523  std::string newExpr = __newExpr.str().c_str();
524 
526  boost::for_each( endExpr, [&isInserted, &geofileString, &newExpr]( boost::regex _expr)
527  {
528  boost::match_results<std::string::iterator> itExpr;
529  bool exprIsFound = boost::regex_search(geofileString.begin(), geofileString.end(), itExpr, _expr);
530  if (exprIsFound)
531  {
532  geofileString.insert(itExpr[0].first, newExpr.begin(), newExpr.end() );
533  isInserted = true;
534  return;
535  }
536  });
537  if (!isInserted)
538  geofileString.insert(geofileString.end(), newExpr.begin(), newExpr.end() );
539  }
540  }
541 
543  std::ostringstream __newAlgo, __newList, __newBgm;
544 
545  if (aniso)
546  __newAlgo << "Field[" << i << "] = MinAniso;\n";
547  else
548  __newAlgo << "Field[" << i << "] = Min;\n";
549 
550  __newList << "Field[" << i << "].FieldsList = {";
551 
552  size_t j=1;
553  for (; j< posfiles.size(); j++)
554  __newList << j << ", ";
555  __newList << j <<"};\n";
556 
557  if (posfiles.size() > 1)
558  __newBgm << "Background Field = " << i << ";\n";
559  else
560  __newBgm << "Background Field = " << i-1 << ";\n";
561 
562  std::string newAlgo = __newAlgo.str().c_str();
563  std::string newList = __newList.str().c_str();
564  std::string newBgm = __newBgm.str().c_str();
565 
566  std::list<std::string> newEndExpr;
567  if (posfiles.size() > 1)
568  {
569  newEndExpr.push_back(newAlgo);
570  newEndExpr.push_back(newList);
571  newEndExpr.push_back(newBgm);
572  }
573  else
574  newEndExpr.push_back(newBgm);
575 
576  std::list<std::string>::iterator itNew = newEndExpr.begin();
577 
580  isInserted = false;
581  for (std::list<boost::regex>::iterator itExpr = endExpr.begin(); itExpr != endExpr.end(); itExpr++)
582  {
583  boost::match_results<std::string::iterator> itPlace;
584  bool exprIsFound = boost::regex_search(geofileString.begin(), geofileString.end(), itPlace, *itExpr);
585  if (exprIsFound)
586  geofileString = boost::regex_replace(geofileString, *itExpr, *itNew);
587  else
588  {
589  std::list<boost::regex>::iterator itNextExpr = itExpr;
590  itNextExpr++;
591  for ( ; itNextExpr != endExpr.end(); itNextExpr++)
592  {
593  bool nextIsFound = boost::regex_search(geofileString.begin(), geofileString.end(), itPlace, *itNextExpr);
594  if (nextIsFound)
595  {
596  geofileString.insert(itPlace[0].first, (*itNew).begin(), (*itNew).end() );
597  isInserted = true;
598  break;
599  }
600  }
601  if (!isInserted)
602  geofileString.insert(geofileString.end(), (*itNew).begin(), (*itNew).end() );
603  }
604  itNew++;
605  }
606 
608  std::ofstream newGeofile;
609 
610  // Find name of the geofile (without extension), from the complete path
611  size_t geoNameBegin = geofile.find_last_of("/");
612  size_t extensionBegin = geofile.find(".geo") - 1;
613  std::string geofileName = geofile.substr( geoNameBegin+1, extensionBegin - geoNameBegin);
614 
615  if( boost::filesystem::exists((boost::format( "./new-%1%.geo" ) % geofileName ).str()) )
616  remove((boost::format( "./new-%1%.geo" ) % geofileName ).str().c_str() );
617 
618  newGeofile.open( (boost::format( "./new-%1%.geo" ) % geofileName ).str().c_str() );
619  newGeofile << geofileString;
620  newGeofile.close();
621 
622  std::string newAccessGeofile = (boost::format( "./new-%1%.geo" ) % geofileName ).str();
623 
624  return newAccessGeofile;
625  }
626 
627  template<int Dim,
628  int Order,
629  int OrderGeo>
630  std::string
631  MeshAdaptation<Dim,
632  Order,
633  OrderGeo>::buildAdaptedMesh(std::string geofile, std::string name, std::vector<std::string> posfiles,
634  bool aniso)
635  {
636  // Find name of the geofile (without extension), from the complete path
637  size_t geoNameBegin = geofile.find_last_of("/");
638  size_t extensionBegin = geofile.find(".geo") - 1;
639  std::string geofileName = geofile.substr( geoNameBegin+1, extensionBegin - geoNameBegin);
640 
641  std::string prefix = (boost::format( "./%1%" ) % geofileName ).str();
642  std::string mshFormat = "msh";
643  std::string newMeshName = (boost::format( "%1%%2%.%3%" ) % prefix % name % mshFormat ).str();
644 
645  int nbPosfiles = posfiles.size();
646 
647  if ( !mpi::environment::initialized() || ( mpi::environment::initialized() && Environment::worldComm().globalRank() == Environment::worldComm().masterRank() ) )
648  {
650 #if FEELPP_HAS_GMSH
651 #if defined(FEELPP_HAS_GMSH_H)
652 
653  std::string geofileNameWE = (boost::format( "%1%.geo" ) % geofileName).str();
654  std::string exeName = "";
655 
656  // Load geofile (.geo) and post processing (.pos) files
657  int argcGmsh = nbPosfiles + 2;
658  char **argvGmsh = new char * [argcGmsh]; // geofile + posfiles
659 
661  char *argExe = new char[exeName.size() + 1];
662  copy(exeName.begin(), exeName.end(), argExe);
663  argExe[exeName.size()] = '\0';
664  argvGmsh[0] = argExe;
665 
666 
668  char *argGeo = new char[geofileNameWE.size() + 1];
669  copy(geofileNameWE.begin(), geofileNameWE.end(), argGeo);
670  argGeo[geofileNameWE.size()] = '\0';
671  argvGmsh[1] = argGeo;
672 
674  for (int i=0; i<nbPosfiles; i++)
675  {
676  int j = i + 2; // Shift (0 = exe, 1 = geo)
677  char *argPos = new char[posfiles[i].size() + 1];
678  copy(posfiles[i].begin(), posfiles[i].end(), argPos);
679  argPos[posfiles[i].size()] = '\0';
680  argvGmsh[j] = argPos;
681  }
682 
684  GmshInitialize(argcGmsh, argvGmsh);
685 
686 
687  GmshSetOption("Mesh", "Algorithm", 5.);
688  ::GModel *newGmshModel = new GModel();
689 
690  //CTX::instance()->mesh.secondOrderIncomplete = 0;
691  //CTX::instance()->mesh.secondOrderLinear = 1; // has to 1 to work
692 
693  // Retreive list of files
694  for (unsigned int i = 0; i < CTX::instance()->files.size(); i++)
695  {
696  LOG(INFO) << "[MeshAdaptation] loaded files : " << CTX::instance()->files[i] << "\n";
697  MergeFile(CTX::instance()->files[i]);
698  }
699 
700  // /* Create Fields from PView list */
701  ::FieldManager* myFieldManager = newGmshModel->getFields();
702  std::list<int> idList;
703 
704  for (unsigned int i = 0; i < ::PView::list.size(); i++)
705  {
706  ::PView *v = ::PView::list[i];
707  if (v->getData()->hasModel(::GModel::current()))
708  {
709  Msg::Error("Cannot use view based on current mesh for background mesh: you might"
710  " want to use a list-based view (.pos file) instead");
711  return 0;
712  }
713 
714  LOG(INFO) << "[MeshAdaptation] PostView : " << v->getData()->getFileName() << "\n";
715 
717  int id = myFieldManager->newId();
718  myFieldManager->newField(id, "PostView");
719 
720  ::Field *f = myFieldManager->get(id);
721  f->options["IView"]->numericalValue(i);
722  idList.push_back(id);
723  }
724 
725  // /* Create minAniso field, with all the posfiles (intersection) */
726  int id = myFieldManager->newId();
727 
728  if (aniso)
729  myFieldManager->newField(id, "MinAniso");
730  else
731  myFieldManager->newField(id, "Min");
732  ::Field *f = myFieldManager->get(id);
733 
735  assert( f->options["FieldsList"]->list().size() == 0);
737  //std::copy(idList.begin(), idList.end(), std::back_inserter(f->options["FieldsList"]->list()) );
738  f->options["FieldsList"]->list(idList);
739 
740  // /* Now create the adapted mesh */
741 
743  myFieldManager->setBackgroundFieldId(id);
744  f = myFieldManager->get(myFieldManager->getBackgroundField());
745 
747  CTX::instance()->mesh.algo2d = ALGO_2D_BAMG;
748  if (Dim == 3)
749  CTX::instance()->mesh.algo3d = ALGO_3D_MMG3D;
750 
751  newGmshModel->deleteMesh(); //Delete current mesh
752  newGmshModel->mesh(Dim);
753 
754  LOG(INFO) << "[MeshAdaptation] New mesh built : " << newMeshName << "\n";
755  newGmshModel->writeMSH(newMeshName);
756  LOG(INFO) << "[MeshAdaptation] vertices : " << newGmshModel->getNumMeshVertices() << "\n";
757  LOG(INFO) << "[MeshAdaptation] elements : " << newGmshModel->getNumMeshElements() << "\n";
758 
760  CTX::instance()->files.erase(CTX::instance()->files.begin(), CTX::instance()->files.end());
762  ::PView::list.erase(::PView::list.begin(), ::PView::list.end() );
763 
764  // Cleanup memory
765  newGmshModel->destroy();
766  delete newGmshModel;
767  delete[] argGeo;
768  delete[] argExe;
769  delete[] argvGmsh;
770  GmshFinalize();
771 
772 #else
773  std::string newAccessGeofile = createAdaptedGeo(geofile, name, posfiles, aniso);
774 
775  // Execute gmsh command to generate new mesh from new geofile built
776  std::ostringstream __str;
777  __str << "gmsh" << " " << newAccessGeofile << " "
778  << "-o " << newMeshName << " ";
779  for (int p=0; p<nbPosfiles; p++)
780  {
781  __str << posfiles[p] << " ";
782  }
783  if (aniso)
784  {
785  __str << "-algo bamg" << " ";
786  if( Dim == 3)
787  __str << "-algo mmg3d" << " ";
788  }
789  __str << "-" << Dim;
790 
791  ::system(__str.str().c_str());
792 #endif
793 #else
794  throw std::invalid_argument("Gmsh is not available on this system");
795 #endif
796  }
797 
798  if ( mpi::environment::initialized() )
799  {
800  Environment::worldComm().barrier();
801  LOG(INFO) << "Broadcast adapted mesh file " << newMeshName << " to all other mpi processes\n";
802  mpi::broadcast( Environment::worldComm().globalComm(), newMeshName, 0 );
803  }
804 
805  return newMeshName;
806  }
807 
808  template<int Dim,
809  int Order,
810  int OrderGeo>
811  void
812  MeshAdaptation<Dim,
813  Order,
814  OrderGeo>::computeMetric(const double tol, const double hMin, const double hMax,
815  const matrixN_type & hessianMatrix, matrixN_type & metric, double & maxEigenvalue, int dofId)
816  {
817  using namespace Feel::vf;
818 
819  Eigen::EigenSolver< matrixN_type > eigenSolver;
820  eigenSolver.compute(hessianMatrix);
821  vectorN_type eigenvalues = eigenSolver.eigenvalues().array().abs();
822 
823  /***** Use of Frey formula - bamg user manual
824  eigenvalues *= 1./(tol*norm);
825  ******/
826 
827  /***** Use of Alauzet formula - p35 Metric-Based Anisotropic Mesh Adaptation ****/
828  //eigenvalues *= 1./(tol*norm);
829  eigenvalues *= 1./tol;
830  eigenvalues *= 0.5*(Dim/(double) (Dim+1))*(Dim/(double) (Dim+1));
831  for (int i=0; i<Dim; i++)
832  {
833  eigenvalues(i)=std::min(std::max(eigenvalues(i),1/(hMax*hMax)), 1/(hMin*hMin));
834  }
835  /******/
836 
837  matrixN_type S = matrixN_type::Zero();
838  for (int i=0; i<Dim; i++)
839  S(i,i)=eigenvalues(i); //1/math::sqrt(eigenvalues(i));
840 
841  matrixN_type R = matrixN_type::Zero();
842  for (int j=0; j<Dim; j++)
843  for (int i=0; i<Dim; i++)
844  R(i,j)=real((eigenSolver.eigenvectors())(i,j)); // because eigenvectors are complex
845 
846  metric = (R * S) * R.transpose();
847  maxEigenvalue = eigenvalues.maxCoeff();
848 
849  for(int i=0; i<Dim; i++)
850  M_measures[dofId](i) = 1.0/math::sqrt(eigenvalues(i));
851 
852  M_directions[dofId] = R;
853  }
854 
855  template<int Dim,
856  int Order,
857  int OrderGeo>
858  std::string
859  MeshAdaptation<Dim,
860  Order,
861  OrderGeo>::adaptMeshHess1(element_type& var, const mesh_ptrtype& mesh, double hMin, double hMax, std::string name,
862  std::string geofile, double tol, bool aniso)
863  {
864  using namespace Feel::vf;
865 
866  p0_space_ptrtype P0h = p0_space_type::New( mesh ); //P0 space
867  p1_space_ptrtype P1h = p1_space_type::New( mesh ); //P1 space
868 
869  // Initialize M_measures and M_directions size
870  M_measures.resize(P1h->nLocalDof());
871  M_directions.resize(P1h->nLocalDof());
872 
873  // Store measure on each point
874  int bbItemSize;
875  if (aniso)
876  bbItemSize = Dim*Dim;
877  else
878  bbItemSize = 1;
879 
880  //cout << "n_bb_item=" << bbItemSize << endl;
881  std::vector<p1_element_type> bbNewMap(bbItemSize, P1h->element()); // map of adapted measures
882 
883  // // Store max and min of solution U
884  // double max_val = U.max();
885  // double min_val = U.min();
886  // double cutoff = 1.e-5;
887  // double norm = std::max(std::max(fabs(max_val),fabs(min_val)), cutoff);
888 
889  // ******* From U in P1 space -> approximate components of hessian matrix on P1 space ******** //
890  // Project u on P1 space
891  auto varP1 = P1h->element();
892  varP1 = vf::project( P1h, elements(mesh), idv(var) );
893 
894  // First derivatives of U (P1) are P0
895  std::vector<p0_element_type> dvard1P0(Dim, P0h->element());
896  std::vector<p1_element_type> dvard1P1(Dim, P1h->element());
897  for (int i=0; i<Dim; i++)
898  {
899  dvard1P0[i] = vf::project(P0h, elements(mesh), gradv(varP1)(0,i) );
900  dvard1P1[i] = element_div( sum(P1h, idv(dvard1P0[i])*meas()), sum(P1h, meas() ) ); //L2 proj -> P1
901  }
902 
903  // Second derivatives of proj_P1(U) are P0
904  std::vector<p0_element_type> hessP0(Dim*Dim, P0h->element());
905  std::vector<p1_element_type> hessP1(Dim*Dim, P1h->element());
906  for (int i=0; i<Dim; i++)
907  {
908  for (int j=0; j<Dim; j++)
909  {
910  // Hessian components are entered column by column
911  hessP0[i+j*Dim] = vf::project(P0h, elements(mesh), gradv( dvard1P1[j] )(0,i) );
912  hessP1[i+j*Dim] = element_div( sum(P1h, idv(hessP0[i+j*Dim])*meas()), sum(P1h, meas() ) ); //L2 proj -> P1
913  }
914  }
915  // ******************************************************************************************* //
916 
917  auto dofptItP1 = P1h->dof()->dofPointBegin();
918  auto dofptEnP1 = P1h->dof()->dofPointEnd();
919 
920  // double hsizeMin = (1.0/1000)*meshSize;
921  // std::cout << "hsize min = " << hsizeMin << std::endl;
922  // double hsizeMax = 1000*meshSize;
923  // std::cout << "hsize max = " << hsizeMax << std::endl;
924 
925  for ( ; dofptItP1 != dofptEnP1; dofptItP1++)
926  {
927  auto dofptCoordP1 = boost::get<0>(*dofptItP1); //dofptItP1->get<0>();
928  auto dofptIdP1 = boost::get<1>(*dofptItP1); //dofptItP1->get<1>();
929 
930  matrixN_type hessianMatrix;
931  Eigen::EigenSolver< matrixN_type > eigenSolver;
932 
933  // Associate each dof with the P1 hessian matrix projection
934  for (int i=0; i<Dim; i++)
935  {
936  for (int j=0; j<Dim; j++)
937  {
938  hessianMatrix(i,j) = (hessP1[i+j*Dim])[dofptIdP1];
939  }
940  }
941 
942  double maxEigenvalue;
943  matrixN_type metrics;
944  computeMetric(tol, hMin, hMax, hessianMatrix, metrics, maxEigenvalue, dofptIdP1);
945 
946  if (aniso)
947  {
948  for (int j=0; j<Dim; j++)
949  {
950  for (int i=0; i<Dim; i++)
951  {
952  (bbNewMap[i+j*Dim])[dofptIdP1] = metrics(i,j);
953  }
954  }
955  }
956  else
957  {
958  auto newHsize = 1.0/math::sqrt( maxEigenvalue);
959  (bbNewMap[0])[dofptIdP1] = newHsize;
960  }
961  }
962 
964  std::string posfileName;
965  if (aniso)
966  posfileName = createPosfileAnisotropic(name, bbNewMap, mesh);
967  else
968  posfileName = createPosfile(name, bbNewMap[0], mesh);
969 
970  return posfileName;
971 
972  }
973 
974  template<int Dim,
975  int Order,
976  int OrderGeo>
977  std::string
978  MeshAdaptation<Dim,
979  Order,
980  OrderGeo>::adaptMeshHess2(element_type& var, const mesh_ptrtype& mesh, double hMin, double hMax,
981  std::string name, std::string geofile, double tol, bool aniso)
982  {
983  return adaptMeshHess2(var, mesh, hMin, hMax, name, geofile, tol, aniso, Feel::mpl::bool_< isP1 >() );
984  }
985 
986  template<int Dim,
987  int Order,
988  int OrderGeo>
989  std::string
990  MeshAdaptation<Dim,
991  Order,
992  OrderGeo>::adaptMeshHess2(element_type& var, const mesh_ptrtype& mesh, double hMin, double hMax,
993  std::string name, std::string geofile, double tol, bool aniso,
994  mpl::bool_<true>)
995  {
996  return adaptMeshHess1(var, mesh, hMin, hMax, name, geofile, tol, aniso);
997  }
998 
999  template<int Dim,
1000  int Order,
1001  int OrderGeo>
1002  std::string
1003  MeshAdaptation<Dim,
1004  Order,
1005  OrderGeo>::adaptMeshHess2(element_type& var, const mesh_ptrtype& mesh, double hMin, double hMax,
1006  std::string name, std::string geofile, double tol, bool aniso,
1007  mpl::bool_<false>)
1008  {
1009  using namespace Feel::vf;
1010 
1011  p1_space_ptrtype P1h = p1_space_type::New( mesh ); //P1 space
1012  p1vec_space_ptrtype P1hvec = p1vec_space_type::New( mesh ); //P1 (vectorial) space
1013 
1014  // Initialize M_measures and M_directions size
1015  M_measures.resize(P1h->nLocalDof());
1016  M_directions.resize(P1h->nLocalDof());
1017 
1018 
1019  // Define P(k-2) space from Order parameter
1020  typedef bases<Lagrange<Order-2,Scalar, Discontinuous> > p_km2_basis_type;
1021  typedef FunctionSpace<mesh_type, p_km2_basis_type> p_km2_space_type;
1022  typedef boost::shared_ptr<p_km2_space_type> p_km2_space_ptrtype;
1023  typedef typename p_km2_space_type::element_type p_km2_element_type;
1024 
1025  p_km2_space_ptrtype Pkm2H = p_km2_space_type::New( mesh );
1026 
1027  // Store measure on each point
1028  int bbItemSize; // = (aniso) ? Dim*Dim : 1;
1029  if (aniso)
1030  bbItemSize = Dim*Dim;
1031  else
1032  bbItemSize = 1;
1033 
1034  //cout << "n_bb_item=" << bbItemSize << endl;
1035  std::vector<p1_element_type> bbNewMap(bbItemSize, P1h->element()); // map of adapted measures
1036 
1037  // Store max and min of solution U
1038  // double max_val = U.max();
1039  // double min_val = U.min();
1040  // double cutoff = 1.e-5;
1041  // double norm = std::max(std::max(fabs(max_val),fabs(min_val)), cutoff);
1042 
1043  // Second derivatives of U are P(k-2)
1044  std::vector<p_km2_element_type> dvard2Pkm2(Dim*Dim, Pkm2H->element() );
1045  for (int i=0; i<Dim; i++)
1046  {
1047  for (int j=0; j<Dim; j++)
1048  {
1049  dvard2Pkm2[i+j*Dim] = vf::project( Pkm2H, elements(mesh),hessv(var)(i,j));
1050  }
1051  }
1052 
1053  p1_element_type V(P1h, "V");
1054 
1057  std::vector<p1_element_type> dvard2P1(Dim*Dim, P1h->element());
1058  //auto l2proj = opProjection( _domainSpace=P1h, _imageSpace=P1h );
1059 
1060  for(int i=0; i<Dim; i++)
1061  {
1062  for (int j=0; j<Dim; j++)
1063  {
1064  //dvard2P1[i+j*Dim] = l2proj->project( _expr=idv( dvard2Pkm2[i+j*Dim]) );
1065  dvard2P1[i+j*Dim] = vf::project( P1h, elements(mesh), idv( dvard2Pkm2[i+j*Dim]) );
1066  }
1067  }
1068  // ******************************************************************************************* //
1069 
1070  auto dofptItP1 = P1h->dof()->dofPointBegin();
1071  auto dofptEnP1 = P1h->dof()->dofPointEnd();
1072 
1073  // double hsizeMin = (1.0/10000)*meshSize;
1074  // std::cout << "hsize min = " << hsizeMin << std::endl;
1075  // double hsizeMax = 10000*meshSize;
1076  // std::cout << "hsize max = " << hsizeMax << std::endl;
1077 
1078  std::vector<p1_element_type> hessianComponents(4);
1079 
1080  for ( ; dofptItP1 != dofptEnP1; dofptItP1++)
1081  {
1082  auto dofptCoordP1 = boost::get<0>(*dofptItP1); //dofptItP1->get<0>();
1083  auto dofptIdP1 = boost::get<1>(*dofptItP1); //dofptItP1->get<1>();
1084 
1085  matrixN_type hessianMatrix;
1086 
1087  // Associate each dof with the P1 hessian matrix projection
1088  for (int i=0; i<Dim; i++)
1089  {
1090  for (int j=0; j<Dim; j++)
1091  {
1092  hessianMatrix(i,j) = (dvard2P1[i+j*Dim])[dofptIdP1];
1093  }
1094  }
1095 
1096  double maxEigenvalue;
1097  matrixN_type metrics;
1098  computeMetric(tol, hMin, hMax, hessianMatrix, metrics, maxEigenvalue, dofptIdP1);
1099 
1100  if (aniso)
1101  {
1102  for (int j=0; j<Dim; j++)
1103  {
1104  for (int i=0; i<Dim; i++)
1105  {
1106  (bbNewMap[i+j*Dim])[dofptIdP1] = metrics(i,j);
1107  }
1108  }
1109  }
1110  else
1111  {
1112  auto newHsize = 1.0/math::sqrt( maxEigenvalue);
1113  (bbNewMap[0])[dofptIdP1] = newHsize;
1114 
1115  }
1116  }
1117 
1119  std::string posfileName;
1120  if (aniso)
1121  posfileName = createPosfileAnisotropic(name, bbNewMap, mesh);
1122  else
1123  posfileName = createPosfile(name, bbNewMap[0], mesh);
1124 
1125  return posfileName;
1126 
1127  }
1128 
1129 }
1130 #endif // __MESHADAPTATION_HPP

Generated on Fri Oct 25 2013 14:24:19 for Feel++ by doxygen 1.8.4