Select Git revision
Options.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Christophe Geuzaine authoredOptions in the GUI are now applied instantenuously. Please let me know if you see anything weird in the GUI.
Christophe Geuzaine authoredOptions in the GUI are now applied instantenuously. Please let me know if you see anything weird in the GUI.
meshMetric.cpp 32.77 KiB
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
#include "meshMetric.h"
#include "meshGFaceOptimize.h"
#include "Context.h"
#include "GEntity.h"
#include "GModel.h"
#include "gmshLevelset.h"
#include "MElementOctree.h"
#include "OS.h"
#include <algorithm>
meshMetric::meshMetric(GModel *gm)
{
hasAnalyticalMetric = false;
_dim = gm->getDim();
std::map<MElement *, MElement *> newP;
std::map<MElement *, MElement *> newD;
if(_dim == 2) {
for(GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); ++fit) {
for(std::size_t i = 0; i < (*fit)->getNumMeshElements(); i++) {
MElement *e = (*fit)->getMeshElement(i);
MElement *copy = e->copy(_vertexMap, newP, newD);
_elements.push_back(copy);
}
}
}
else if(_dim == 3) {
for(GModel::riter rit = gm->firstRegion(); rit != gm->lastRegion(); ++rit) {
for(std::size_t i = 0; i < (*rit)->getNumMeshElements(); i++) {
MElement *e = (*rit)->getMeshElement(i);
MElement *copy = e->copy(_vertexMap, newP, newD);
_elements.push_back(copy);
}
}
}
_octree = new MElementOctree(_elements);
buildVertexToElement(_elements, _adj);
}
meshMetric::meshMetric(std::vector<MElement *> elements)
{
hasAnalyticalMetric = false;
_dim = elements[0]->getDim();
std::map<MElement *, MElement *> newP;
std::map<MElement *, MElement *> newD;
for(std::size_t i = 0; i < elements.size(); i++) {
MElement *e = elements[i];
MElement *copy = e->copy(_vertexMap, newP, newD);
_elements.push_back(copy);
}
_octree = new MElementOctree(_elements);
buildVertexToElement(_elements, _adj);
}
void meshMetric::addMetric(int technique, simpleFunction<double> *fct,
const std::vector<double> ¶meters)
{
needMetricUpdate = true;
int metricNumber = setOfMetrics.size();
setOfFcts[metricNumber] = fct;
setOfParameters[metricNumber] = parameters; setOfTechniques[metricNumber] = technique;
if(fct->hasDerivatives()) hasAnalyticalMetric = true;
computeMetric(metricNumber);
}
void meshMetric::updateMetrics()
{
if(!setOfMetrics.size()) {
Msg::Error("Can't intersect metrics, no metric registered");
return;
}
v2t_cont ::iterator it = _adj.begin();
for(; it != _adj.end(); it++) {
MVertex *ver = it->first;
_nodalMetrics[ver] = setOfMetrics[0][ver];
_nodalSizes[ver] = setOfSizes[0][ver];
if(setOfMetrics.size() > 1)
for(std::size_t i = 1; i < setOfMetrics.size(); i++) {
_nodalMetrics[ver] =
(_dim == 3) ? intersection_conserve_mostaniso(_nodalMetrics[ver],
setOfMetrics[i][ver]) :
intersection_conserve_mostaniso_2d(
_nodalMetrics[ver], setOfMetrics[i][ver]);
_nodalSizes[ver] = std::min(_nodalSizes[ver], setOfSizes[i][ver]);
}
}
needMetricUpdate = false;
}
void meshMetric::exportInfo(const char *fileendname)
{
if(needMetricUpdate) updateMetrics();
std::stringstream sg, sm, sl, sh, shm;
sg << "meshmetric_gradients_" << fileendname;
sm << "meshmetric_metric_" << fileendname;
sl << "meshmetric_levelset_" << fileendname;
sh << "meshmetric_hessian_" << fileendname;
shm << "meshmetric_hessianmat_" << fileendname;
std::ofstream out_grad(sg.str().c_str());
std::ofstream out_metric(sm.str().c_str());
std::ofstream out_ls(sl.str().c_str());
std::ofstream out_hess(sh.str().c_str());
std::ofstream out_hessmat(shm.str().c_str());
out_grad << "View \"ls_gradient\"{" << std::endl;
out_metric << "View \"metric\"{" << std::endl;
out_ls << "View \"ls\"{" << std::endl;
out_hess << "View \"hessian\"{" << std::endl;
out_hessmat << "View \"hessian_mat\"{" << std::endl;
std::vector<MElement *>::iterator itelem = _elements.begin();
std::vector<MElement *>::iterator itelemen = _elements.end();
for(; itelem != itelemen; itelem++) {
MElement *e = *itelem;
if(e->getDim() == 2) {
out_metric << "TT(";
out_grad << "VT(";
out_ls << "ST(";
out_hess << "ST(";
out_hessmat << "TT(";
}
else {
out_metric << "TS(";
out_grad << "VS(";
out_ls << "SS(";
out_hess << "SS(";
out_hessmat << "TS(";
} for(std::size_t i = 0; i < e->getNumVertices(); i++) {
MVertex *ver = e->getVertex(i);
out_metric << ver->x() << "," << ver->y() << "," << ver->z();
out_grad << ver->x() << "," << ver->y() << "," << ver->z();
out_ls << ver->x() << "," << ver->y() << "," << ver->z();
out_hess << ver->x() << "," << ver->y() << "," << ver->z();
out_hessmat << ver->x() << "," << ver->y() << "," << ver->z();
if(i != e->getNumVertices() - 1) {
out_metric << ",";
out_grad << ",";
out_ls << ",";
out_hess << ",";
out_hessmat << ",";
}
else {
out_metric << "){";
out_grad << "){";
out_ls << "){";
out_hess << "){";
out_hessmat << "){";
}
}
for(std::size_t i = 0; i < e->getNumVertices(); i++) {
MVertex *ver = e->getVertex(i);
out_ls << vals[ver];
out_hess << (hessians[ver](0, 0) + hessians[ver](1, 1) +
hessians[ver](2, 2));
if(i == (e->getNumVertices() - 1)) {
out_ls << "};" << std::endl;
out_hess << "};" << std::endl;
}
else {
out_ls << ",";
out_hess << ",";
}
for(int k = 0; k < 3; k++) {
out_grad << grads[ver](k);
if((k == 2) && (i == (e->getNumVertices() - 1)))
out_grad << "};" << std::endl;
else
out_grad << ",";
for(int l = 0; l < 3; l++) {
out_metric << _nodalMetrics[ver](k, l);
out_hessmat << hessians[ver](k, l);
if((k == 2) && (l == 2) && (i == (e->getNumVertices() - 1))) {
out_metric << "};" << std::endl;
out_hessmat << "};" << std::endl;
}
else {
out_metric << ",";
out_hessmat << ",";
}
}
}
}
}
out_grad << "};" << std::endl;
out_metric << "};" << std::endl;
out_ls << "};" << std::endl;
out_hess << "};" << std::endl;
out_hessmat << "};" << std::endl;
out_grad.close();
out_metric.close();
out_ls.close();
out_hess.close();
out_hessmat.close();
}
meshMetric::~meshMetric()
{ if(_octree) delete _octree;
for(std::size_t i = 0; i < _elements.size(); i++) delete _elements[i];
}
void meshMetric::computeValues()
{
v2t_cont ::iterator it = _adj.begin();
while(it != _adj.end()) {
std::vector<MElement *> lt = it->second;
MVertex *ver = it->first;
vals[ver] = (*_fct)(ver->x(), ver->y(), ver->z());
it++;
}
}
// Determines set of vertices to use for least squares
std::vector<MVertex *> getLSBlob(std::size_t minNbPt, v2t_cont::iterator it,
v2t_cont &adj)
{
// static const double RADFACT = 3;
//
// SPoint3 p0 = it->first->point(); // Vertex coordinates (center of circle)
//
// double rad = 0.;
// for (std::vector<MElement*>::iterator itEl = it->second.begin(); itEl !=
// it->second.end(); itEl++)
// rad = std::max(rad,(*itEl)->getOuterRadius());
// rad *= RADFACT;
std::vector<MVertex *> vv(1, it->first),
bvv = vv; // Vector of vertices in blob and in boundary of blob
do {
std::set<MVertex *> nbvv; // Set of vertices in new boundary
for(std::vector<MVertex *>::iterator itBV = bvv.begin(); itBV != bvv.end();
itBV++) { // For each boundary vertex...
std::vector<MElement *> &adjBV = adj[*itBV];
for(std::vector<MElement *>::iterator itBVEl = adjBV.begin();
itBVEl != adjBV.end(); itBVEl++) {
for(std::size_t iV = 0; iV < (*itBVEl)->getNumVertices();
iV++) { // ... look for adjacent vertices...
MVertex *v = (*itBVEl)->getVertex(iV);
// if ((find(vv.begin(),vv.end(),v) == vv.end()) &&
// (p0.distance(v->point()) <= rad)) nbvv.insert(v);
if(find(vv.begin(), vv.end(), v) == vv.end())
nbvv.insert(v); // ... and add them in the new boundary if they are
// not already in the blob
}
}
}
if(nbvv.empty())
bvv.clear();
else {
bvv.assign(nbvv.begin(), nbvv.end());
vv.insert(vv.end(), nbvv.begin(), nbvv.end());
}
// } while (!bvv.empty());
} while(vv.size() < minNbPt); // Repeat until min. number of points is reached
return vv;
}
// Compute derivatives and second order derivatives using least squares
// 2D LS system: a_i0*x^2+a_i1*x*y+a_i2*y^2+a_i3*x+a_i4*y+a_i5=b_i
// 3D LS system:
// a_i0*x^2+a_i1*x*y+a_i2*x*z+a_i3*y^2+a_i4*y*z+a_i5*z^2+a_i6*x+a_i7*y+a_i8*z+a_i9=b_i
void meshMetric::computeHessian()
{
std::size_t sysDim = (_dim == 2) ? 6 : 10;
std::size_t minNbPtBlob = 3 * sysDim;
for(v2t_cont::iterator it = _adj.begin(); it != _adj.end(); it++) {
MVertex *ver = it->first;
std::vector<MVertex *> vv = getLSBlob(minNbPtBlob, it, _adj);
fullMatrix<double> A(vv.size(), sysDim), ATA(sysDim, sysDim);
fullVector<double> b(vv.size()), ATb(sysDim), coeffs(sysDim);
for(std::size_t i = 0; i < vv.size(); i++) {
const double &x = vv[i]->x(), &y = vv[i]->y(), &z = vv[i]->z();
if(_dim == 2) {
A(i, 0) = x * x;
A(i, 1) = x * y;
A(i, 2) = y * y;
A(i, 3) = x;
A(i, 4) = y;
A(i, 5) = 1.;
}
else {
A(i, 0) = x * x;
A(i, 1) = x * y;
A(i, 2) = x * z;
A(i, 3) = y * y;
A(i, 4) = y * z;
A(i, 5) = z * z;
A(i, 6) = x;
A(i, 7) = y;
A(i, 8) = z;
A(i, 9) = 1.;
}
b(i) = vals[vv[i]];
}
ATA.gemm(A, A, 1., 0., true, false);
A.multWithATranspose(b, 1., 0., ATb);
ATA.luSolve(ATb, coeffs);
const double &x = ver->x(), &y = ver->y(), &z = ver->z();
double d2udx2, d2udy2, d2udz2, d2udxy, d2udxz, d2udyz, dudx, dudy, dudz;
if(_dim == 2) {
d2udx2 = 2. * coeffs(0);
d2udy2 = 2. * coeffs(2);
d2udz2 = 0.;
d2udxy = coeffs(1);
d2udxz = 0.;
d2udyz = 0.;
dudx = d2udx2 * x + d2udxy * y + coeffs(3);
dudy = d2udxy * x + d2udy2 * y + coeffs(4);
dudz = 0.;
}
else {
d2udx2 = 2. * coeffs(0);
d2udy2 = 2. * coeffs(3);
d2udz2 = 2. * coeffs(5);
d2udxy = coeffs(1);
d2udxz = coeffs(2);
d2udyz = coeffs(4);
dudx = d2udx2 * x + d2udxy * y + d2udxz * z + coeffs(6);
dudy = d2udxy * x + d2udy2 * y + d2udyz * z + coeffs(7);
dudz = d2udxz * x + d2udyz * y + d2udz2 * z + coeffs(8);
}
double duNorm = sqrt(dudx * dudx + dudy * dudy + dudz * dudz);
if(duNorm == 0. || _technique == meshMetric::HESSIAN ||
_technique == meshMetric::EIGENDIRECTIONS ||
_technique == meshMetric::EIGENDIRECTIONS_LINEARINTERP_H)
duNorm = 1.;
grads[ver] = SVector3(dudx / duNorm, dudy / duNorm, dudz / duNorm);
hessians[ver](0, 0) = d2udx2;
hessians[ver](0, 1) = d2udxy;
hessians[ver](0, 2) = d2udxz;
hessians[ver](1, 0) = d2udxy;
hessians[ver](1, 1) = d2udy2;
hessians[ver](1, 2) = d2udyz;
hessians[ver](2, 0) = d2udxz;
hessians[ver](2, 1) = d2udyz; hessians[ver](2, 2) = d2udz2;
}
}
void meshMetric::computeMetricLevelSet(MVertex *ver, SMetric3 &hessian,
SMetric3 &metric, double &size, double x,
double y, double z)
{
double signed_dist;
SVector3 gr;
if(ver) {
signed_dist = vals[ver];
gr = grads[ver];
hessian = hessians[ver];
}
else {
signed_dist = (*_fct)(x, y, z);
_fct->gradient(x, y, z, gr(0), gr(1), gr(2));
_fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
hessian(2, 1), hessian(2, 2));
}
double dist = fabs(signed_dist);
SMetric3 H;
double norm = gr(0) * gr(0) + gr(1) * gr(1) + gr(2) * gr(2);
if(dist < _e && norm != 0.0) {
double h = hmin * (hmax / hmin - 1) * dist / _e + hmin;
double C = 1. / (h * h) - 1. / (hmax * hmax);
H(0, 0) += C * gr(0) * gr(0) / norm;
H(1, 1) += C * gr(1) * gr(1) / norm;
H(2, 2) += C * gr(2) * gr(2) / norm;
H(1, 0) = H(0, 1) = C * gr(1) * gr(0) / norm;
H(2, 0) = H(0, 2) = C * gr(2) * gr(0) / norm;
H(2, 1) = H(1, 2) = C * gr(2) * gr(1) / norm;
}
fullMatrix<double> V(3, 3);
fullVector<double> S(3);
H.eig(V, S);
double lambda1, lambda2, lambda3;
lambda1 = S(0);
lambda2 = S(1);
lambda3 = (_dim == 3) ? S(2) : 1.;
SVector3 t1(V(0, 0), V(1, 0), V(2, 0));
SVector3 t2(V(0, 1), V(1, 1), V(2, 1));
SVector3 t3(V(0, 2), V(1, 2), V(2, 2));
size =
std::min(std::min(1 / sqrt(lambda1), 1 / sqrt(lambda2)), 1 / sqrt(lambda3));
metric = SMetric3(lambda1, lambda2, lambda3, t1, t2, t3);
}
void meshMetric::computeMetricHessian(MVertex *ver, SMetric3 &hessian,
SMetric3 &metric, double &size, double x,
double y, double z)
{
SVector3 gr;
if(ver != NULL) {
gr = grads[ver];
hessian = hessians[ver];
}
else if(ver == NULL) {
_fct->gradient(x, y, z, gr(0), gr(1), gr(2));
_fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
hessian(2, 1), hessian(2, 2)); }
double _epsilonP = 1.;
double hminP = 1.e-12;
double hmaxP = 1.e+12;
fullMatrix<double> V(3, 3);
fullVector<double> S(3);
hessian.eig(V, S);
double lambda1 =
std::min(std::max(fabs(S(0)) / _epsilonP, 1. / (hmaxP * hmaxP)),
1. / (hminP * hminP));
double lambda2 =
std::min(std::max(fabs(S(1)) / _epsilonP, 1. / (hmaxP * hmaxP)),
1. / (hminP * hminP));
double lambda3 =
(_dim == 3) ?
std::min(std::max(fabs(S(2)) / _epsilonP, 1. / (hmaxP * hmaxP)),
1. / (hminP * hminP)) :
1.;
SVector3 t1(V(0, 0), V(1, 0), V(2, 0));
SVector3 t2(V(0, 1), V(1, 1), V(2, 1));
SVector3 t3 =
(_dim == 3) ? SVector3(V(0, 2), V(1, 2), V(2, 2)) : SVector3(0., 0., 1.);
size =
std::min(std::min(1 / sqrt(lambda1), 1 / sqrt(lambda2)), 1 / sqrt(lambda3));
metric = SMetric3(lambda1, lambda2, lambda3, t1, t2, t3);
}
void meshMetric::computeMetricFrey(MVertex *ver, SMetric3 &hessian,
SMetric3 &metric, double &size, double x,
double y, double z)
{
double signed_dist;
SVector3 gr;
if(ver) {
signed_dist = vals[ver];
gr = grads[ver];
hessian = hessians[ver];
}
else {
signed_dist = (*_fct)(x, y, z);
_fct->gradient(x, y, z, gr(0), gr(1), gr(2));
_fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
hessian(2, 1), hessian(2, 2));
}
double dist = fabs(signed_dist);
SMetric3 H(1. / (hmax * hmax));
double norm = gr(0) * gr(0) + gr(1) * gr(1) + gr(2) * gr(2);
if(dist < _e && norm != 0.0) {
double h = hmin * (hmax / hmin - 1.0) * dist / _e + hmin;
double C = 1. / (h * h) - 1. / (hmax * hmax);
double kappa = hessian(0, 0) + hessian(1, 1) + hessian(2, 2);
double epsGeom = 4.0 * 3.14 * 3.14 / (kappa * _np * _np);
H(0, 0) += C * gr(0) * gr(0) / (norm) + hessian(0, 0) / epsGeom;
H(1, 1) += C * gr(1) * gr(1) / (norm) + hessian(1, 1) / epsGeom;
H(2, 2) += C * gr(2) * gr(2) / (norm) + hessian(2, 2) / epsGeom;
H(1, 0) = H(0, 1) = C * gr(1) * gr(0) / (norm) + hessian(1, 0) / epsGeom;
H(2, 0) = H(0, 2) = C * gr(2) * gr(0) / (norm) + hessian(2, 0) / epsGeom;
H(2, 1) = H(1, 2) = C * gr(2) * gr(1) / (norm) + hessian(2, 1) / epsGeom;
}
fullMatrix<double> V(3, 3);
fullVector<double> S(3); H.eig(V, S);
double lambda1, lambda2, lambda3;
lambda1 = S(0);
lambda2 = S(1);
lambda3 = (_dim == 3) ? S(2) : 1.;
if(dist < _e) {
lambda1 = std::min(std::max(fabs(S(0)) / _epsilon, 1. / (hmax * hmax)),
1. / (hmin * hmin));
lambda2 = std::min(std::max(fabs(S(1)) / _epsilon, 1. / (hmax * hmax)),
1. / (hmin * hmin));
lambda3 = (_dim == 3) ?
std::min(std::max(fabs(S(2)) / _epsilon, 1. / (hmax * hmax)),
1. / (hmin * hmin)) :
1.;
}
SVector3 t1(V(0, 0), V(1, 0), V(2, 0));
SVector3 t2(V(0, 1), V(1, 1), V(2, 1));
SVector3 t3(V(0, 2), V(1, 2), V(2, 2));
size =
std::min(std::min(1 / sqrt(lambda1), 1 / sqrt(lambda2)), 1 / sqrt(lambda3));
metric = SMetric3(lambda1, lambda2, lambda3, t1, t2, t3);
}
void meshMetric::computeMetricEigenDir(MVertex *ver, SMetric3 &hessian,
SMetric3 &metric, double &size, double x,
double y, double z)
{
double signed_dist;
SVector3 gVec;
if(ver) {
signed_dist = vals[ver];
gVec = grads[ver];
hessian = hessians[ver];
}
else {
signed_dist = (*_fct)(x, y, z);
_fct->gradient(x, y, z, gVec(0), gVec(1), gVec(2));
_fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
hessian(2, 1), hessian(2, 2));
}
double dist = fabs(signed_dist);
const double metric_value_hmax = 1. / (hmax * hmax);
const double gMag = gVec.norm(), invGMag = 1. / gMag;
if(signed_dist < _e && signed_dist > _e_moins && gMag != 0.0) {
const double metric_value_hmin = 1. / (hmin * hmin);
const SVector3 nVec = invGMag * gVec; // Unit normal vector
double lambda_n =
0.; // Eigenvalues of metric for normal & tangential directions
if(_technique == meshMetric::EIGENDIRECTIONS_LINEARINTERP_H) {
// const double h_dist = hmin + ((hmax-hmin)/_e)*dist;
// // Characteristic element size in the normal direction - linear interp
// between hmin and hmax lambda_n = 1./(h_dist*h_dist);
double beta = CTX::instance()->mesh.smoothRatio;
double h_dista = std::min((hmin + (dist * log(beta))), hmax);
lambda_n = 1. / (h_dista * h_dista);
}
else if(_technique == meshMetric::EIGENDIRECTIONS) {
const double maximum_distance =
(signed_dist > 0.) ? _e : fabs(_e_moins); // ... or linear interpolation
// between 1/h_min^2 and
// 1/h_max^2
lambda_n = metric_value_hmin +
((metric_value_hmax - metric_value_hmin) / maximum_distance) * dist;
}
std::vector<SVector3> tVec; // Unit tangential vectors
std::vector<double> kappa; // Curvatures
if(_dim == 2) { // 2D curvature formula: cf. R. Goldman, "Curvature formulas
// for implicit curves and surfaces", Computer Aided
// Geometric Design 22 (2005), pp. 632–658
kappa.resize(2);
kappa[0] =
fabs(-gVec(1) * (-gVec(1) * hessian(0, 0) + gVec(0) * hessian(0, 1)) +
gVec(0) * (-gVec(1) * hessian(1, 0) + gVec(0) * hessian(1, 1))) *
pow(invGMag, 3);
kappa[1] = 1.;
tVec.resize(2);
tVec[0] = SVector3(-nVec(1), nVec(0), 0.);
tVec[1] = SVector3(0., 0., 1.);
}
else { // 3D curvature formula: cf. A.G. Belyaev, A.A. Pasko and T.L. Kunii,
// "Ridges and Ravines on Implicit Surfaces," CGI, pp.530-535,
// Computer Graphics International 1998 (CGI'98), 1998
fullMatrix<double> ImGG(3, 3);
ImGG(0, 0) = 1. - gVec(0) * gVec(0);
ImGG(0, 1) = -gVec(0) * gVec(1);
ImGG(0, 2) = -gVec(0) * gVec(2);
ImGG(1, 0) = -gVec(1) * gVec(0);
ImGG(1, 1) = 1. - gVec(1) * gVec(1);
ImGG(1, 2) = -gVec(1) * gVec(2);
ImGG(2, 0) = -gVec(2) * gVec(0);
ImGG(2, 1) = -gVec(2) * gVec(1);
ImGG(2, 2) = 1. - gVec(2) * gVec(2);
fullMatrix<double> hess(3, 3);
hessian.getMat(hess);
fullMatrix<double> gN(3, 3); // Gradient of unit normal
gN.gemm(ImGG, hess, 1., 0.);
gN.scale(invGMag);
fullMatrix<double> eigVecL(3, 3), eigVecR(3, 3);
fullVector<double> eigValRe(3), eigValIm(3);
gN.eig(eigValRe, eigValIm, eigVecL, eigVecR,
false); // Eigendecomp. of gradient of unit normal
kappa.resize(3); // Store abs. val. of eigenvalues (= principal curvatures
// only in non-normal directions)
kappa[0] = fabs(eigValRe(0));
kappa[1] = fabs(eigValRe(1));
kappa[2] = fabs(eigValRe(2));
tVec.resize(3); // Store normalized eigenvectors (= principal directions
// only in non-normal directions)
tVec[0] = SVector3(eigVecR(0, 0), eigVecR(1, 0), eigVecR(2, 0));
tVec[0].normalize();
tVec[1] = SVector3(eigVecR(0, 1), eigVecR(1, 1), eigVecR(2, 1));
tVec[1].normalize();
tVec[2] = SVector3(eigVecR(0, 2), eigVecR(1, 2), eigVecR(2, 2));
tVec[2].normalize();
std::vector<double> tVecDotNVec(3); // Store dot products with normal
// vector to look for normal direction
tVecDotNVec[0] = fabs(dot(tVec[0], nVec));
tVecDotNVec[1] = fabs(dot(tVec[1], nVec));
tVecDotNVec[2] = fabs(dot(tVec[2], nVec));
const int i_N = max_element(tVecDotNVec.begin(), tVecDotNVec.end()) -
tVecDotNVec.begin(); // Index of normal dir. = max. dot
// products (close to 0. in
// tangential dir.)
kappa.erase(kappa.begin() + i_N); // Remove normal dir.
tVec.erase(tVec.begin() + i_N);
}
const double invh_t0 = (_np * kappa[0]) / 6.283185307,
invh_t1 = (_np * kappa[1]) /
6.283185307; // Inverse of tangential size 0
const double lambda_t0 = invh_t0 * invh_t0, lambda_t1 = invh_t1 * invh_t1;
const double lambdaP_n = std::min(std::max(lambda_n, metric_value_hmax),
metric_value_hmin); // Truncate eigenvalues
const double lambdaP_t0 =
std::min(std::max(lambda_t0, metric_value_hmax), metric_value_hmin);
const double lambdaP_t1 =
(_dim == 2) ?
1. :
std::min(std::max(lambda_t1, metric_value_hmax), metric_value_hmin);
metric =
SMetric3(lambdaP_n, lambdaP_t0, lambdaP_t1, nVec, tVec[0], tVec[1]);
const double h_n = 1. / sqrt(lambdaP_n), h_t0 = 1. / sqrt(lambdaP_t0),
h_t1 = 1. / sqrt(lambdaP_t1);
size = std::min(std::min(h_n, h_t0), h_t1);
}
else { // isotropic metric !
SMetric3 mymetric(metric_value_hmax);
metric = mymetric;
size = hmax;
}
}
void meshMetric::computeMetricIsoLinInterp(MVertex *ver, SMetric3 &hessian,
SMetric3 &metric, double &size,
double x, double y, double z)
{
double signed_dist;
SVector3 gr;
if(ver) {
signed_dist = vals[ver];
gr = grads[ver];
hessian = hessians[ver];
}
else {
signed_dist = (*_fct)(x, y, z);
_fct->gradient(x, y, z, gr(0), gr(1), gr(2));
_fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
hessian(2, 1), hessian(2, 2));
}
double dist = fabs(signed_dist);
double norm = gr.normalize();
size = hmax; // the characteristic element size in all directions - linear
// interp between hmin and hmax
if(norm != 0.) {
if((signed_dist >= 0) && (signed_dist < _e))
size = hmin + ((hmax - hmin) / _e) * dist;
else if((signed_dist < 0) && (signed_dist > _e_moins))
size = hmin - ((hmax - hmin) / _e_moins) * dist;
}
double lambda = 1. / size / size;
metric(0, 0) = lambda;
metric(0, 1) = 0.;
metric(0, 2) = 0.;
metric(1, 0) = 0.;
metric(1, 1) = lambda;
metric(1, 2) = 0.;
metric(2, 0) = 0.;
metric(2, 1) = 0.;
metric(2, 2) = (_dim == 3) ? lambda : 1.;
}
// this function scales the mesh metric in order
// to reach a target number of elements
// We know that the number of elements in the final
// mesh will be (assuming M_e the metric at centroid of element e)
// N = \sum_e \sqrt {\det (M_e)} V_e
// where V_e is the volume of e
// assuming that N_{target} = K N, we have// K N = K \sum_e \sqrt {\det (M_e)} V_e
// = \sum_e \sqrt {K^2 \det (M_e)} V_e
// = \sum_e \sqrt {\det (K^{2/d} M_e)} V_e
// where d is the dimension of the problem.
// This means that the metric should be scaled by K^{2/d} where
// K is N_target / N
void meshMetric::scaleMetric(int nbElementsTarget, nodalMetricTensor &nmt)
{
// compute N
double N = 0;
for(std::size_t i = 0; i < _elements.size(); i++) {
MElement *e = _elements[i];
SMetric3 m1 = nmt[e->getVertex(0)];
SMetric3 m2 = nmt[e->getVertex(1)];
SMetric3 m3 = nmt[e->getVertex(2)];
if(_dim == 2) {
SMetric3 m = interpolation(m1, m2, m3, 0.3333, 0.3333);
N += sqrt(m.determinant()) * e->getVolume() * 4. / sqrt(3.0); // 3.0
}
else {
SMetric3 m4 = nmt[e->getVertex(3)];
SMetric3 m = interpolation(m1, m2, m3, m4, 0.25, 0.25, 0.25);
N += sqrt(m.determinant()) * e->getVolume() * 12. / sqrt(2.0); // 4.0;
}
}
double scale = pow((double)nbElementsTarget / N, 2.0 / _dim);
for(nodalMetricTensor::iterator it = nmt.begin(); it != nmt.end(); ++it) {
if(_dim == 3) {
it->second *= scale;
}
else {
it->second(0, 0) *= scale;
it->second(1, 0) *= scale;
it->second(1, 1) *= scale;
}
SMetric3 &m = it->second;
fullMatrix<double> V(3, 3);
fullVector<double> S(3);
m.eig(V, S);
S(0) = std::min(std::max(S(0), 1 / (hmax * hmax)), 1 / (hmin * hmin));
S(1) = std::min(std::max(S(1), 1 / (hmax * hmax)), 1 / (hmin * hmin));
if(_dim == 3)
S(2) = std::min(std::max(S(2), 1 / (hmax * hmax)), 1 / (hmin * hmin));
SVector3 t1(V(0, 0), V(1, 0), V(2, 0));
SVector3 t2(V(0, 1), V(1, 1), V(2, 1));
SVector3 t3(V(0, 2), V(1, 2), V(2, 2));
m = SMetric3(S(0), S(1), S(2), t1, t2, t3);
}
}
void meshMetric::computeMetric(int metricNumber)
{
_fct = setOfFcts[metricNumber];
std::vector<double> parameters = setOfParameters[metricNumber];
int technique = setOfTechniques[metricNumber];
hmin = parameters.size() >= 3 ? parameters[1] : CTX::instance()->mesh.lcMin;
hmax = parameters.size() >= 3 ? parameters[2] : CTX::instance()->mesh.lcMax;
_e = parameters[0];
_e_moins = (parameters.size() >= 5) ? parameters[4] : -parameters[0];
if(_e_moins > 0.) _e_moins *= -1.;
_epsilon = parameters[0];
_technique = (MetricComputationTechnique)technique;
_np = (parameters.size() >= 4) ? parameters[3] : 15.;
computeValues();
computeHessian();
for(v2t_cont::iterator it = _adj.begin(); it != _adj.end(); it++) {
MVertex *ver = it->first; SMetric3 hessian, metric;
double size;
switch(_technique) {
case(LEVELSET): computeMetricLevelSet(ver, hessian, metric, size); break;
case(HESSIAN): computeMetricHessian(ver, hessian, metric, size); break;
case(FREY): computeMetricFrey(ver, hessian, metric, size); break;
case(EIGENDIRECTIONS):
computeMetricEigenDir(ver, hessian, metric, size);
break;
case(EIGENDIRECTIONS_LINEARINTERP_H):
computeMetricEigenDir(ver, hessian, metric, size);
break;
case(ISOTROPIC_LINEARINTERP_H):
computeMetricIsoLinInterp(ver, hessian, metric, size);
break;
}
setOfSizes[metricNumber].insert(std::make_pair(ver, size));
setOfMetrics[metricNumber].insert(std::make_pair(ver, metric));
}
if(_technique == HESSIAN) scaleMetric(_epsilon, setOfMetrics[metricNumber]);
}
double meshMetric::operator()(double x, double y, double z, GEntity *ge)
{
if(needMetricUpdate) updateMetrics();
if(!setOfMetrics.size()) {
std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
throw;
}
SPoint3 xyz(x, y, z), uvw;
double initialTol = MElement::getTolerance();
MElement::setTolerance(1.e-4);
MElement *e = _octree->find(x, y, z, _dim);
MElement::setTolerance(initialTol);
double value = 0.;
if(e) {
e->xyz2uvw(xyz, uvw);
double *val = new double[e->getNumVertices()];
for(std::size_t i = 0; i < e->getNumVertices(); i++) {
val[i] = _nodalSizes[e->getVertex(i)];
}
value = e->interpolate(val, uvw[0], uvw[1], uvw[2]);
delete[] val;
}
else {
Msg::Warning("point %g %g %g not found, looking for nearest node", x, y, z);
double minDist = 1.e100;
for(nodalField::iterator it = _nodalSizes.begin(); it != _nodalSizes.end();
it++) {
const double dist = xyz.distance(it->first->point());
if(dist <= minDist) {
minDist = dist;
value = it->second;
}
}
}
return value;
}
void meshMetric::operator()(double x, double y, double z, SMetric3 &metr,
GEntity *ge)
{
if(needMetricUpdate) {
updateMetrics();
}
if(!setOfMetrics.size()) {
std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
throw; }
metr = SMetric3(1.e-22);
// RECOMPUTE MESH METRIC AT XYZ
if(hasAnalyticalMetric) {
int nbMetrics = setOfMetrics.size();
std::vector<SMetric3> newSetOfMetrics(nbMetrics);
for(int iMetric = 0; iMetric < nbMetrics; iMetric++) {
_fct = setOfFcts[iMetric];
_technique = (MetricComputationTechnique)setOfTechniques[iMetric];
if(_fct->hasDerivatives()) {
SMetric3 hessian, metric;
double size;
switch(_technique) {
case(LEVELSET):
computeMetricLevelSet(NULL, hessian, metric, size, x, y, z);
break;
case(HESSIAN):
computeMetricHessian(NULL, hessian, metric, size, x, y, z);
break;
case(FREY):
computeMetricFrey(NULL, hessian, metric, size, x, y, z);
break;
case(EIGENDIRECTIONS):
computeMetricEigenDir(NULL, hessian, metric, size, x, y, z);
break;
case(EIGENDIRECTIONS_LINEARINTERP_H):
computeMetricEigenDir(NULL, hessian, metric, size, x, y, z);
break;
case(ISOTROPIC_LINEARINTERP_H):
computeMetricIsoLinInterp(NULL, hessian, metric, size, x, y, z);
break;
}
newSetOfMetrics[iMetric] = metric;
}
else {
// find other metrics here
SMetric3 metric;
SPoint3 xyz(x, y, z), uvw;
double initialTol = MElement::getTolerance();
MElement::setTolerance(1.e-4);
MElement *e = _octree->find(x, y, z, _dim);
MElement::setTolerance(initialTol);
if(e) {
e->xyz2uvw(xyz, uvw);
SMetric3 m1 = setOfMetrics[iMetric][e->getVertex(0)];
SMetric3 m2 = setOfMetrics[iMetric][e->getVertex(1)];
SMetric3 m3 = setOfMetrics[iMetric][e->getVertex(2)];
if(_dim == 2)
metric = interpolation(m1, m2, m3, uvw[0], uvw[1]);
else {
SMetric3 m4 = setOfMetrics[iMetric][e->getVertex(3)];
metric = interpolation(m1, m2, m3, m4, uvw[0], uvw[1], uvw[2]);
}
newSetOfMetrics[iMetric] = metric;
}
else {
Msg::Warning("point %g %g %g not found, looking for nearest node", x,
y, z);
}
}
}
// intersect metrics here
metr = newSetOfMetrics[0];
for(int i = 1; i < nbMetrics; i++)
metr = intersection_conserve_mostaniso(metr, newSetOfMetrics[i]);
}
// INTERPOLATE DISCRETE MESH METRIC
else {
SPoint3 xyz(x, y, z), uvw; double initialTol = MElement::getTolerance();
MElement::setTolerance(1.e-4);
MElement *e = _octree->find(x, y, z, _dim);
MElement::setTolerance(initialTol);
if(e) {
e->xyz2uvw(xyz, uvw);
SMetric3 m1 = _nodalMetrics[e->getVertex(0)];
SMetric3 m2 = _nodalMetrics[e->getVertex(1)];
SMetric3 m3 = _nodalMetrics[e->getVertex(2)];
if(_dim == 2)
metr = interpolation(m1, m2, m3, uvw[0], uvw[1]);
else {
SMetric3 m4 = _nodalMetrics[e->getVertex(3)];
metr = interpolation(m1, m2, m3, m4, uvw[0], uvw[1], uvw[2]);
}
}
else {
Msg::Warning("point %g %g %g not found, looking for nearest node", x, y,
z);
double minDist = 1.e100;
for(nodalMetricTensor::iterator it = _nodalMetrics.begin();
it != _nodalMetrics.end(); it++) {
const double dist = xyz.distance(it->first->point());
if(dist <= minDist) {
minDist = dist;
metr = it->second;
}
}
}
}
}
double meshMetric::getLaplacian(MVertex *v)
{
MVertex *vNew = _vertexMap[v->getNum()];
std::map<MVertex *, SMetric3>::const_iterator it = hessians.find(vNew);
SMetric3 h = it->second;
return h(0, 0) + h(1, 1) + h(2, 2);
}
SVector3 meshMetric::getGradient(MVertex *v)
{
MVertex *vNew = _vertexMap[v->getNum()];
std::map<MVertex *, SVector3>::const_iterator it = grads.find(vNew);
SVector3 gr = it->second;
return gr;
}