30 #ifndef __MESHADAPTATION_HPP
31 #define __MESHADAPTATION_HPP 1
37 #include <feel/feeldiscr/functionspace.hpp>
45 #include <feel/feelcore/feel.hpp>
47 #include <feel/feeldiscr/projector.hpp>
49 #include <boost/assert.hpp>
50 #include <boost/regex.hpp>
54 #if defined( FEELPP_HAS_GMSH_H )
60 #include <StringUtils.h>
63 #include <PViewData.h>
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)
84 static const bool isP1 = (Order == 1);
87 typedef Eigen::Matrix<double, Dim, 1> vectorN_type;
88 typedef Eigen::Matrix<double, Dim, Dim> matrixN_type;
91 typedef Backend<double> backend_type;
93 typedef boost::shared_ptr<backend_type> backend_ptrtype;
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;
104 typedef Simplex<Dim,OrderGeo> convex_type;
106 typedef Mesh<convex_type> mesh_type;
108 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
111 typedef Exporter<mesh_type> export_type;
113 typedef boost::shared_ptr<export_type> export_ptrtype;
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;
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;
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;
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;
140 typedef typename mesh_type::point_type point_type;
146 element_type initVar;
147 std::pair<element_type, std::string> initPairVar = std::make_pair( initVar,
"defaultVar");
148 M_defaultVar.push_back(initPairVar);
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);
159 std::vector<vectorN_type>
const measures()
165 std::vector<matrixN_type>
const directions()
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);
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);
181 std::string createAdaptedGeo(std::string geofile, std::string name, std::vector<std::string> posfiles,
bool aniso);
184 std::string buildAdaptedMesh(std::string geofile, std::string name, std::vector<std::string> posfiles ,
bool aniso);
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);
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);
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,
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,
214 BOOST_PARAMETER_MEMBER_FUNCTION(
219 (initMesh,(mesh_ptrtype))
220 (geofile, (std::string))
221 (adaptType, (std::string)))
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))
230 mesh_ptrtype newMesh = this-> adaptMeshImpl( initMesh, geofile, adaptType, var, metric, hMin, hMax, tol);
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;
238 std::vector<vectorN_type> M_measures;
239 std::vector<matrixN_type> M_directions;
247 MeshAdaptation<Dim, Order, OrderGeo>::isP1;
252 typename MeshAdaptation<Dim, Order, OrderGeo>::mesh_ptrtype
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)
260 std::vector<std::string> posfiles;
261 std::string finalName =
"";
262 if ( !(var.front().second ==
"defaultVar") )
264 boost::for_each( var, [&initMesh, &hMin, &hMax, &geofile, &tol, &adaptType, &posfiles, &finalName,
this]
265 (std::pair<element_type, std::string> _var)
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;
274 else if ( !(metric.front().second ==
"defaultMetric") )
276 boost::for_each(metric, [&initMesh, &adaptType, &posfiles, &finalName,
this]
277 (std::pair<std::vector<p1_element_type>, std::string> _metric)
279 std::string newPosfile;
280 if (adaptType ==
"isotropic")
282 BOOST_ASSERT(_metric.first.size()==1);
283 newPosfile = this->createPosfile(_metric.second, _metric.first[0], initMesh);
285 else if( adaptType ==
"anisotropic")
287 BOOST_ASSERT(_metric.first.size()==Dim*Dim);
288 newPosfile = this->createPosfileAnisotropic(_metric.second, _metric.first, initMesh);
290 posfiles.push_back( newPosfile );
291 finalName +=
"-" + _metric.second;
295 std::string newMeshName = this->buildAdaptedMesh( geofile, finalName, posfiles, adaptType==
"anisotropic");
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";
314 OrderGeo>::createPosfile(std::string nameVar,
const p1_element_type& bbNewMap,
const mesh_ptrtype& mesh)
316 using namespace Feel::vf;
318 std::string posFormat =
"pos";
319 std::string mshFormat =
"msh";
321 p1_space_ptrtype P1h = p1_space_type::New( mesh );
323 std::ofstream newPosFile( (boost::format(
"%1%.%2%" ) % nameVar %posFormat ).str() );
325 newPosFile <<
"View \" background mesh \" { \n";
326 auto eltIt = mesh->beginElement();
327 auto eltEnd = mesh->endElement();
328 for ( ; eltIt != eltEnd; eltIt++)
330 std::vector<point_type> eltPoints( eltIt->nPoints() );
335 for (
int i=0; i < eltIt->nPoints(); i++)
337 eltPoints[i] = eltIt->point(i);
338 for (
int j=0; j< Dim; j++)
340 newPosFile << eltPoints[i](j);
341 if ( Dim == 2 || (Dim == 3 && j!=Dim-1) )
349 if (i!= eltIt->nPoints() - 1)
355 for (
size_t k=0; k<eltPoints.size(); k++)
357 auto dofIndex = boost::get<0>( P1h->dof()->localToGlobal( eltIt->id(), k) );
358 newPosFile << bbNewMap[ dofIndex ];
359 if ( k!= eltPoints.size() - 1)
363 newPosFile <<
"};\n";
370 return (boost::format(
"%1%.%2%" ) % nameVar %posFormat ).str();
379 OrderGeo>::createPosfileAnisotropic(std::string nameVar,
const std::vector<p1_element_type>& bbNewMap,
380 const mesh_ptrtype& mesh)
382 using namespace Feel::vf;
384 std::string posFormat =
"pos";
385 std::string mshFormat =
"msh";
387 p1_space_ptrtype P1h = p1_space_type::New( mesh );
389 std::ofstream newPosFile( (boost::format(
"%1%.%2%" ) % nameVar %posFormat ).str() );
391 newPosFile <<
"View \" background mesh \" { \n";
392 auto eltIt = mesh->beginElement();
393 auto eltEnd = mesh->endElement();
394 for ( ; eltIt != eltEnd; eltIt++)
396 std::vector<point_type> eltPoints( eltIt->nPoints() );
404 for (
int i=0; i < eltIt->nPoints(); i++)
406 eltPoints[i] = eltIt->point(i);
407 for (
int j=0; j< Dim; j++)
409 newPosFile << eltPoints[i](j);
410 if ( Dim == 2 || (Dim == 3 && j!=Dim-1) )
415 if (i!= eltIt->nPoints() - 1)
421 for (
size_t k=0; k<eltPoints.size(); k++)
423 auto dofIndex = boost::get<0>( P1h->dof()->localToGlobal( eltIt->id(), k) );
426 newPosFile << (bbNewMap[num++])[ dofIndex ];
429 newPosFile << (bbNewMap[num++])[ dofIndex ] <<
", ";
432 newPosFile << (bbNewMap[num++])[ dofIndex ] <<
", ";
436 newPosFile << (bbNewMap[num++])[ dofIndex ] <<
", ";
437 newPosFile << (bbNewMap[num++])[ dofIndex ] <<
", ";
441 newPosFile << (bbNewMap[num++])[ dofIndex ] <<
", ";
442 newPosFile << (bbNewMap[num++])[ dofIndex ] <<
", ";
443 newPosFile << (bbNewMap[num++])[ dofIndex ] <<
", ";
444 newPosFile << (bbNewMap[num++])[ dofIndex ];
454 if (k!= eltPoints.size() - 1)
457 newPosFile <<
"};\n";
464 return (boost::format(
"%1%.%2%" ) % nameVar %posFormat ).str();
473 OrderGeo>::createAdaptedGeo(std::string geofile, std::string name, std::vector<std::string> posfiles,
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>());
487 int nbPosfiles = posfiles.size();
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:]]*);";
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);
500 bool isInserted =
false;
504 for (; i <= nbPosfiles ; i++)
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 <<
";" ;
510 boost::regex expression1( __expr1.str().c_str() );
511 boost::regex expression2( __expr2.str().c_str() );
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);
518 if (!(fieldIsFoundL1 && fieldIsFoundL2))
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();
526 boost::for_each( endExpr, [&isInserted, &geofileString, &newExpr]( boost::regex _expr)
528 boost::match_results<std::string::iterator> itExpr;
529 bool exprIsFound = boost::regex_search(geofileString.begin(), geofileString.end(), itExpr, _expr);
532 geofileString.insert(itExpr[0].first, newExpr.begin(), newExpr.end() );
538 geofileString.insert(geofileString.end(), newExpr.begin(), newExpr.end() );
543 std::ostringstream __newAlgo, __newList, __newBgm;
546 __newAlgo <<
"Field[" << i <<
"] = MinAniso;\n";
548 __newAlgo <<
"Field[" << i <<
"] = Min;\n";
550 __newList <<
"Field[" << i <<
"].FieldsList = {";
553 for (; j< posfiles.size(); j++)
554 __newList << j <<
", ";
555 __newList << j <<
"};\n";
557 if (posfiles.size() > 1)
558 __newBgm <<
"Background Field = " << i <<
";\n";
560 __newBgm <<
"Background Field = " << i-1 <<
";\n";
562 std::string newAlgo = __newAlgo.str().c_str();
563 std::string newList = __newList.str().c_str();
564 std::string newBgm = __newBgm.str().c_str();
566 std::list<std::string> newEndExpr;
567 if (posfiles.size() > 1)
569 newEndExpr.push_back(newAlgo);
570 newEndExpr.push_back(newList);
571 newEndExpr.push_back(newBgm);
574 newEndExpr.push_back(newBgm);
576 std::list<std::string>::iterator itNew = newEndExpr.begin();
581 for (std::list<boost::regex>::iterator itExpr = endExpr.begin(); itExpr != endExpr.end(); itExpr++)
583 boost::match_results<std::string::iterator> itPlace;
584 bool exprIsFound = boost::regex_search(geofileString.begin(), geofileString.end(), itPlace, *itExpr);
586 geofileString = boost::regex_replace(geofileString, *itExpr, *itNew);
589 std::list<boost::regex>::iterator itNextExpr = itExpr;
591 for ( ; itNextExpr != endExpr.end(); itNextExpr++)
593 bool nextIsFound = boost::regex_search(geofileString.begin(), geofileString.end(), itPlace, *itNextExpr);
596 geofileString.insert(itPlace[0].first, (*itNew).begin(), (*itNew).end() );
602 geofileString.insert(geofileString.end(), (*itNew).begin(), (*itNew).end() );
608 std::ofstream newGeofile;
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);
615 if( boost::filesystem::exists((boost::format(
"./new-%1%.geo" ) % geofileName ).str()) )
616 remove((boost::format(
"./new-%1%.geo" ) % geofileName ).str().c_str() );
618 newGeofile.open( (boost::format(
"./new-%1%.geo" ) % geofileName ).str().c_str() );
619 newGeofile << geofileString;
622 std::string newAccessGeofile = (boost::format(
"./new-%1%.geo" ) % geofileName ).str();
624 return newAccessGeofile;
633 OrderGeo>::buildAdaptedMesh(std::string geofile, std::string name, std::vector<std::string> posfiles,
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);
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();
645 int nbPosfiles = posfiles.size();
647 if ( !mpi::environment::initialized() || ( mpi::environment::initialized() && Environment::worldComm().globalRank() == Environment::worldComm().masterRank() ) )
651 #if defined(FEELPP_HAS_GMSH_H)
653 std::string geofileNameWE = (boost::format(
"%1%.geo" ) % geofileName).str();
654 std::string exeName =
"";
657 int argcGmsh = nbPosfiles + 2;
658 char **argvGmsh =
new char * [argcGmsh];
661 char *argExe =
new char[exeName.size() + 1];
662 copy(exeName.begin(), exeName.end(), argExe);
663 argExe[exeName.size()] =
'\0';
664 argvGmsh[0] = argExe;
668 char *argGeo =
new char[geofileNameWE.size() + 1];
669 copy(geofileNameWE.begin(), geofileNameWE.end(), argGeo);
670 argGeo[geofileNameWE.size()] =
'\0';
671 argvGmsh[1] = argGeo;
674 for (
int i=0; i<nbPosfiles; i++)
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;
684 GmshInitialize(argcGmsh, argvGmsh);
687 GmshSetOption(
"Mesh",
"Algorithm", 5.);
688 ::GModel *newGmshModel =
new GModel();
694 for (
unsigned int i = 0; i < CTX::instance()->files.size(); i++)
696 LOG(INFO) <<
"[MeshAdaptation] loaded files : " << CTX::instance()->files[i] <<
"\n";
697 MergeFile(CTX::instance()->files[i]);
701 ::FieldManager* myFieldManager = newGmshModel->getFields();
702 std::list<int> idList;
704 for (
unsigned int i = 0; i < ::PView::list.size(); i++)
706 ::PView *v = ::PView::list[i];
707 if (v->getData()->hasModel(::GModel::current()))
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");
714 LOG(INFO) <<
"[MeshAdaptation] PostView : " << v->getData()->getFileName() <<
"\n";
717 int id = myFieldManager->newId();
718 myFieldManager->newField(
id,
"PostView");
720 ::Field *f = myFieldManager->get(
id);
721 f->options[
"IView"]->numericalValue(i);
722 idList.push_back(
id);
726 int id = myFieldManager->newId();
729 myFieldManager->newField(
id,
"MinAniso");
731 myFieldManager->newField(
id,
"Min");
732 ::Field *f = myFieldManager->get(
id);
735 assert( f->options[
"FieldsList"]->list().size() == 0);
738 f->options[
"FieldsList"]->list(idList);
743 myFieldManager->setBackgroundFieldId(
id);
744 f = myFieldManager->get(myFieldManager->getBackgroundField());
747 CTX::instance()->mesh.algo2d = ALGO_2D_BAMG;
749 CTX::instance()->mesh.algo3d = ALGO_3D_MMG3D;
751 newGmshModel->deleteMesh();
752 newGmshModel->mesh(Dim);
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";
760 CTX::instance()->files.erase(CTX::instance()->files.begin(), CTX::instance()->files.end());
762 ::PView::list.erase(::PView::list.begin(), ::PView::list.end() );
765 newGmshModel->destroy();
773 std::string newAccessGeofile = createAdaptedGeo(geofile, name, posfiles, aniso);
776 std::ostringstream __str;
777 __str <<
"gmsh" <<
" " << newAccessGeofile <<
" "
778 <<
"-o " << newMeshName <<
" ";
779 for (
int p=0; p<nbPosfiles; p++)
781 __str << posfiles[p] <<
" ";
785 __str <<
"-algo bamg" <<
" ";
787 __str <<
"-algo mmg3d" <<
" ";
791 ::system(__str.str().c_str());
794 throw std::invalid_argument(
"Gmsh is not available on this system");
798 if ( mpi::environment::initialized() )
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 );
814 OrderGeo>::computeMetric(
const double tol,
const double hMin,
const double hMax,
815 const matrixN_type & hessianMatrix, matrixN_type & metric,
double & maxEigenvalue,
int dofId)
817 using namespace Feel::vf;
819 Eigen::EigenSolver< matrixN_type > eigenSolver;
820 eigenSolver.compute(hessianMatrix);
821 vectorN_type eigenvalues = eigenSolver.eigenvalues().array().abs();
829 eigenvalues *= 1./tol;
830 eigenvalues *= 0.5*(Dim/(double) (Dim+1))*(Dim/(
double) (Dim+1));
831 for (
int i=0; i<Dim; i++)
833 eigenvalues(i)=std::min(std::max(eigenvalues(i),1/(hMax*hMax)), 1/(hMin*hMin));
837 matrixN_type S = matrixN_type::Zero();
838 for (
int i=0; i<Dim; i++)
839 S(i,i)=eigenvalues(i);
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));
846 metric = (R * S) * R.transpose();
847 maxEigenvalue = eigenvalues.maxCoeff();
849 for(
int i=0; i<Dim; i++)
850 M_measures[dofId](i) = 1.0/math::sqrt(eigenvalues(i));
852 M_directions[dofId] = R;
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)
864 using namespace Feel::vf;
866 p0_space_ptrtype P0h = p0_space_type::New( mesh );
867 p1_space_ptrtype P1h = p1_space_type::New( mesh );
870 M_measures.resize(P1h->nLocalDof());
871 M_directions.resize(P1h->nLocalDof());
876 bbItemSize = Dim*Dim;
881 std::vector<p1_element_type> bbNewMap(bbItemSize, P1h->element());
891 auto varP1 = P1h->element();
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++)
900 dvard1P1[i] =
element_div( sum(P1h, idv(dvard1P0[i])*meas()), sum(P1h, meas() ) );
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++)
908 for (
int j=0; j<Dim; j++)
912 hessP1[i+j*Dim] =
element_div( sum(P1h, idv(hessP0[i+j*Dim])*meas()), sum(P1h, meas() ) );
917 auto dofptItP1 = P1h->dof()->dofPointBegin();
918 auto dofptEnP1 = P1h->dof()->dofPointEnd();
925 for ( ; dofptItP1 != dofptEnP1; dofptItP1++)
927 auto dofptCoordP1 = boost::get<0>(*dofptItP1);
928 auto dofptIdP1 = boost::get<1>(*dofptItP1);
930 matrixN_type hessianMatrix;
931 Eigen::EigenSolver< matrixN_type > eigenSolver;
934 for (
int i=0; i<Dim; i++)
936 for (
int j=0; j<Dim; j++)
938 hessianMatrix(i,j) = (hessP1[i+j*Dim])[dofptIdP1];
942 double maxEigenvalue;
943 matrixN_type metrics;
944 computeMetric(tol, hMin, hMax, hessianMatrix, metrics, maxEigenvalue, dofptIdP1);
948 for (
int j=0; j<Dim; j++)
950 for (
int i=0; i<Dim; i++)
952 (bbNewMap[i+j*Dim])[dofptIdP1] = metrics(i,j);
958 auto newHsize = 1.0/math::sqrt( maxEigenvalue);
959 (bbNewMap[0])[dofptIdP1] = newHsize;
964 std::string posfileName;
966 posfileName = createPosfileAnisotropic(name, bbNewMap, mesh);
968 posfileName = createPosfile(name, bbNewMap[0], mesh);
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)
983 return adaptMeshHess2(var, mesh, hMin, hMax, name, geofile, tol, aniso, Feel::mpl::bool_< isP1 >() );
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,
996 return adaptMeshHess1(var, mesh, hMin, hMax, name, geofile, tol, aniso);
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,
1009 using namespace Feel::vf;
1011 p1_space_ptrtype P1h = p1_space_type::New( mesh );
1012 p1vec_space_ptrtype P1hvec = p1vec_space_type::New( mesh );
1015 M_measures.resize(P1h->nLocalDof());
1016 M_directions.resize(P1h->nLocalDof());
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;
1025 p_km2_space_ptrtype Pkm2H = p_km2_space_type::New( mesh );
1030 bbItemSize = Dim*Dim;
1035 std::vector<p1_element_type> bbNewMap(bbItemSize, P1h->element());
1044 std::vector<p_km2_element_type> dvard2Pkm2(Dim*Dim, Pkm2H->element() );
1045 for (
int i=0; i<Dim; i++)
1047 for (
int j=0; j<Dim; j++)
1053 p1_element_type V(P1h,
"V");
1057 std::vector<p1_element_type> dvard2P1(Dim*Dim, P1h->element());
1060 for(
int i=0; i<Dim; i++)
1062 for (
int j=0; j<Dim; j++)
1070 auto dofptItP1 = P1h->dof()->dofPointBegin();
1071 auto dofptEnP1 = P1h->dof()->dofPointEnd();
1078 std::vector<p1_element_type> hessianComponents(4);
1080 for ( ; dofptItP1 != dofptEnP1; dofptItP1++)
1082 auto dofptCoordP1 = boost::get<0>(*dofptItP1);
1083 auto dofptIdP1 = boost::get<1>(*dofptItP1);
1085 matrixN_type hessianMatrix;
1088 for (
int i=0; i<Dim; i++)
1090 for (
int j=0; j<Dim; j++)
1092 hessianMatrix(i,j) = (dvard2P1[i+j*Dim])[dofptIdP1];
1096 double maxEigenvalue;
1097 matrixN_type metrics;
1098 computeMetric(tol, hMin, hMax, hessianMatrix, metrics, maxEigenvalue, dofptIdP1);
1102 for (
int j=0; j<Dim; j++)
1104 for (
int i=0; i<Dim; i++)
1106 (bbNewMap[i+j*Dim])[dofptIdP1] = metrics(i,j);
1112 auto newHsize = 1.0/math::sqrt( maxEigenvalue);
1113 (bbNewMap[0])[dofptIdP1] = newHsize;
1119 std::string posfileName;
1121 posfileName = createPosfileAnisotropic(name, bbNewMap, mesh);
1123 posfileName = createPosfile(name, bbNewMap[0], mesh);
1130 #endif // __MESHADAPTATION_HPP
elements_type const & elements() const
Definition: elements.hpp:355
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
ElementType element_div(ElementType const &v1, ElementType const &v2)
Definition: functionspace.hpp:6268