Skip to content
Snippets Groups Projects
Commit 4d2b0ea1 authored by Amaury Johnen's avatar Amaury Johnen
Browse files

AnalyseCurvedMesh development (works for triangle and tetra but is not yet... 

AnalyseCurvedMesh development (works for triangle and tetra but is not yet optimized and is still incomplete, edit PluginManager.cpp to activate the plugin)
parent 9829b70d
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Showing
with 2293 additions and 52 deletions
Geo/MElement.h
+ 4
− 0
View file @ 4d2b0ea1
......@@ -15,6 +15,7 @@
#include "MEdge.h"
#include "MFace.h"
#include "polynomialBasis.h"
#include "JacobianBasis.h"
#include "Gauss.h"
class GFace;
......@@ -205,6 +206,9 @@ class MElement
// get the function space for the element
virtual const polynomialBasis* getFunctionSpace(int order=-1) const { return 0; }
// get the function space for the jacobian of the element
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const { return 0; }
// return the interpolating nodal shape functions evaluated at point
// (u,v,w) in parametric coordinates (if order == -1, use the
// polynomial order of the element)
......
Geo/MTetrahedron.cpp
+ 34
− 0
View file @ 4d2b0ea1
......@@ -137,6 +137,40 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
return 0;
}
const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
int nv = getNumVolumeVertices();
if ((nv == 0) && (o == -1)) {
switch (order) {
case 1: return JacobianBases::find(MSH_TET_4);
case 2: return JacobianBases::find(MSH_TET_10);
case 3: return JacobianBases::find(MSH_TET_20);
case 4: return JacobianBases::find(MSH_TET_34);
case 5: return JacobianBases::find(MSH_TET_52);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
else {
switch (order) {
case 1: return JacobianBases::find(MSH_TET_4);
case 2: return JacobianBases::find(MSH_TET_10);
case 3: return JacobianBases::find(MSH_TET_20);
case 4: return JacobianBases::find(MSH_TET_35);
case 5: return JacobianBases::find(MSH_TET_56);
case 6: return JacobianBases::find(MSH_TET_84);
case 7: return JacobianBases::find(MSH_TET_120);
case 8: return JacobianBases::find(MSH_TET_165);
case 9: return JacobianBases::find(MSH_TET_220);
case 10: return JacobianBases::find(MSH_TET_286);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
return 0;
}
int MTetrahedron10::getNumEdgesRep(){ return 6 * CTX::instance()->mesh.numSubEdges; }
int MTetrahedronN::getNumEdgesRep(){ return 6 * CTX::instance()->mesh.numSubEdges; }
......
Geo/MTetrahedron.h
+ 1
− 0
View file @ 4d2b0ea1
......@@ -133,6 +133,7 @@ class MTetrahedron : public MElement {
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
......
Geo/MTriangle.cpp
+ 39
− 0
View file @ 4d2b0ea1
......@@ -104,6 +104,45 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
return 0;
}
const JacobianBasis* MTriangle::getJacobianFuncSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
int nf = getNumFaceVertices();
if ((nf == 0) && (o == -1)) {
switch (order) {
case 1: return JacobianBases::find(MSH_TRI_3);
case 2: return JacobianBases::find(MSH_TRI_6);
case 3: return JacobianBases::find(MSH_TRI_9);
case 4: return JacobianBases::find(MSH_TRI_12);
case 5: return JacobianBases::find(MSH_TRI_15I);
case 6: return JacobianBases::find(MSH_TRI_18);
case 7: return JacobianBases::find(MSH_TRI_21I);
case 8: return JacobianBases::find(MSH_TRI_24);
case 9: return JacobianBases::find(MSH_TRI_27);
case 10: return JacobianBases::find(MSH_TRI_30);
default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
}
}
else {
switch (order) {
case 1: return JacobianBases::find(MSH_TRI_3);
case 2: return JacobianBases::find(MSH_TRI_6);
case 3: return JacobianBases::find(MSH_TRI_10);
case 4: return JacobianBases::find(MSH_TRI_15);
case 5: return JacobianBases::find(MSH_TRI_21);
case 6: return JacobianBases::find(MSH_TRI_28);
case 7: return JacobianBases::find(MSH_TRI_36);
case 8: return JacobianBases::find(MSH_TRI_45);
case 9: return JacobianBases::find(MSH_TRI_55);
case 10: return JacobianBases::find(MSH_TRI_66);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
return 0;
}
int MTriangleN::getNumEdgesRep(){ return 3 * CTX::instance()->mesh.numSubEdges; }
int MTriangle6::getNumEdgesRep(){ return 3 * CTX::instance()->mesh.numSubEdges; }
......
Geo/MTriangle.h
+ 1
− 0
View file @ 4d2b0ea1
......@@ -126,6 +126,7 @@ class MTriangle : public MElement {
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
......
Numeric/CMakeLists.txt
+ 1
− 0
View file @ 4d2b0ea1
......@@ -7,6 +7,7 @@ set(SRC
Numeric.cpp
fullMatrix.cpp
polynomialBasis.cpp
JacobianBasis.cpp
GaussQuadratureLin.cpp
GaussQuadratureTri.cpp
GaussQuadratureQuad.cpp
......
Numeric/JacobianBasis.cpp 0 → 100644
+ 1679
− 0
View file @ 4d2b0ea1
This diff is collapsed.
Numeric/JacobianBasis.h 0 → 100644
+ 237
− 0
View file @ 4d2b0ea1
// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#ifndef _JACOBIAN_BASIS_H_
#define _JACOBIAN_BASIS_H_
#include <math.h>
#include <map>
#include <vector>
#include "fullMatrix.h"
#include <iostream>
class JacobianBasis
{
public:
fullMatrix<double> exposants; //exposants of Bezier FF
fullMatrix<double> points; //sampling point
fullMatrix<double> matrixLag2Bez;
fullMatrix<double> gradShapeMatX;
fullMatrix<double> gradShapeMatY;
fullMatrix<double> gradShapeMatZ;
};
class JacobianBases
{
private:
static std::map<int, JacobianBasis> fs;
public:
static const JacobianBasis *find(int);
};
//// presently those function spaces are only for simplices and quads;
//// should be extended to other elements like hexes
//class polynomialBasis
//{
// public:
// typedef std::vector<std::vector<int> > clCont;
// clCont closures;
// fullMatrix<double> points;
// fullMatrix<double> monomials;
// fullMatrix<double> coefficients;
// int numFaces;
// // for a given face/edge, with both a sign and a rotation,
// // give an ordered list of nodes on this face/edge
// inline const std::vector<int> &getClosure(int id) const // return the closure of dimension dim
// {
// return closures[id];
// }
// inline int getClosureId(int iEl, int iSign=1, int iRot=0) const {
// return iEl + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot;
// }
// inline void evaluateMonomials(double u, double v, double w, double p[]) const
// {
// for (int j = 0; j < monomials.size1(); j++) {
// p[j] = pow_int(u, (int)monomials(j, 0));
// if (monomials.size2() > 1) p[j] *= pow_int(v, (int)monomials(j, 1));
// if (monomials.size2() > 2) p[j] *= pow_int(w, (int)monomials(j, 2));
// }
// }
// inline void f(double u, double v, double w, double *sf) const
// {
// double p[256];
// evaluateMonomials(u, v, w, p);
// for (int i = 0; i < coefficients.size1(); i++) {
// sf[i] = 0;
// for (int j = 0; j < coefficients.size2(); j++) {
// sf[i] += coefficients(i, j) * p[j];
// }
// }
// }
// // I would favour this interface that allows optimizations (not implemented) and is easier to bind
// inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf) {
// double p[256];
// sf.resize (coord.size1(), coefficients.size1());
// for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
// evaluateMonomials(coord(iPoint,0), coord(iPoint,1), coord(iPoint,2), p);
// for (int i = 0; i < coefficients.size1(); i++)
// for (int j = 0; j < coefficients.size2(); j++)
// sf(iPoint,i) += coefficients(i, j) * p[j];
// }
// }
// inline void df(double u, double v, double w, double grads[][3]) const
// {
// switch (monomials.size2()) {
// case 1:
// for (int i = 0; i < coefficients.size1(); i++){
// grads[i][0] = 0;
// grads[i][1] = 0;
// grads[i][2] = 0;
// for(int j = 0; j < coefficients.size2(); j++){
// if (monomials(j, 0) > 0)
// grads[i][0] += coefficients(i, j) *
// pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0);
// }
// }
// break;
// case 2:
// for (int i = 0; i < coefficients.size1(); i++){
// grads[i][0] = 0;
// grads[i][1] = 0;
// grads[i][2] = 0;
// for(int j = 0; j < coefficients.size2(); j++){
// if (monomials(j, 0) > 0)
// grads[i][0] += coefficients(i, j) *
// pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
// pow_int(v, (int)monomials(j, 1));
// if (monomials(j, 1) > 0)
// grads[i][1] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0)) *
// pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1);
// }
// }
// break;
// case 3:
// for (int i = 0; i < coefficients.size1(); i++){
// grads[i][0] = 0;
// grads[i][1] = 0;
// grads[i][2] = 0;
// for(int j = 0; j < coefficients.size2(); j++){
// if (monomials(j, 0) > 0)
// grads[i][0] += coefficients(i, j) *
// pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
// pow_int(v, (int)monomials(j, 1)) *
// pow_int(w, (int)monomials(j, 2));
// if (monomials(j, 1) > 0)
// grads[i][1] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0)) *
// pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1) *
// pow_int(w, (int)monomials(j, 2));
// if (monomials(j, 2) > 0)
// grads[i][2] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0)) *
// pow_int(v, (int)monomials(j, 1)) *
// pow_int(w, (int)(monomials(j, 2) - 1)) * monomials(j, 2);
// }
// }
// break;
// }
// }
// inline void ddf(double u, double v, double w, double hess[][3][3]) const
// {
// switch (monomials.size2()) {
// case 1:
// for (int i = 0; i < coefficients.size1(); i++){
// hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
// hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
// hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
//
// for(int j = 0; j < coefficients.size2(); j++){
// if (monomials(j, 0) > 1) // second derivative !=0
// hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
// monomials(j, 0) * (monomials(j, 0) - 1);
// }
// }
// break;
// case 2:
// for (int i = 0; i < coefficients.size1(); i++){
// hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
// hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
// hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
// for(int j = 0; j < coefficients.size2(); j++){
// if (monomials(j, 0) > 1) // second derivative !=0
// hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
// monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, (int)monomials(j, 1));
// if ((monomials(j, 1) > 0) && (monomials(j, 0) > 0))
// hess[i][0][1] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
// pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
// if (monomials(j, 1) > 1)
// hess[i][1][1] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0)) *
// pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
// }
// hess[i][1][0] = hess[i][0][1];
// }
// break;
// case 3:
// for (int i = 0; i < coefficients.size1(); i++){
// hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
// hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
// hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
// for(int j = 0; j < coefficients.size2(); j++){
// if (monomials(j, 0) > 1)
// hess[i][0][0] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
// pow_int(v, (int)monomials(j, 1)) *
// pow_int(w, (int)monomials(j, 2));
//
// if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0))
// hess[i][0][1] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
// pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
// pow_int(w, (int)monomials(j, 2));
// if ((monomials(j, 0) > 0) && (monomials(j, 2) > 0))
// hess[i][0][2] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
// pow_int(v, (int)monomials(j, 1)) *
// pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
// if (monomials(j, 1) > 1)
// hess[i][1][1] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0)) *
// pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
// pow_int(w, (int)monomials(j, 2));
// if ((monomials(j, 1) > 0) && (monomials(j, 2) > 0))
// hess[i][1][2] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0)) *
// pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
// pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
// if (monomials(j, 2) > 1)
// hess[i][2][2] += coefficients(i, j) *
// pow_int(u, (int)monomials(j, 0)) *
// pow_int(v, (int)monomials(j, 1)) *
// pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
// }
// hess[i][1][0] = hess[i][0][1];
// hess[i][2][0] = hess[i][0][2];
// hess[i][2][1] = hess[i][1][2];
// }
// break;
// }
// }
// static void registerBindings(binding *b);
//};
//
//class polynomialBases
//{
// private:
// static std::map<int, polynomialBasis> fs;
// static std::map<std::pair<int, int>, fullMatrix<double> > injector;
// public :
// static const polynomialBasis *find(int);
// static const fullMatrix<double> &findInjector(int, int);
//};
#endif
Numeric/fullMatrix.cpp
+ 2
− 2
View file @ 4d2b0ea1
......@@ -110,7 +110,7 @@ void fullMatrix<std::complex<double> >::gemm(const fullMatrix<std::complex<doubl
}
template<>
void fullMatrix<double>::mult(const fullVector<double> &x, fullVector<double> &y)
void fullMatrix<double>::mult(const fullVector<double> &x, fullVector<double> &y) const
{
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
double alpha = 1., beta = 0.;
......@@ -120,7 +120,7 @@ void fullMatrix<double>::mult(const fullVector<double> &x, fullVector<double> &y
template<>
void fullMatrix<std::complex<double> >::mult(const fullVector<std::complex<double> > &x,
fullVector<std::complex<double> > &y)
fullVector<std::complex<double> > &y) const
{
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
std::complex<double> alpha = 1., beta = 0.;
......
Numeric/fullMatrix.h
+ 1
− 1
View file @ 4d2b0ea1
......@@ -321,7 +321,7 @@ class fullMatrix
for(int j = 0; j < size2(); j++)
(*this)(i, j) += a*m(i, j);
}
void mult(const fullVector<scalar> &x, fullVector<scalar> &y)
void mult(const fullVector<scalar> &x, fullVector<scalar> &y) const
#if !defined(HAVE_BLAS)
{
y.scale(0.);
......
Plugin/AnalyseCurvedMesh.cpp
+ 284
− 48
View file @ 4d2b0ea1
......@@ -5,14 +5,16 @@
//
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "GmshDefines.h"
#include <stdlib.h>
#include "Gmsh.h"
#include "MTriangle.h"
#include "GmshConfig.h"
#include "GModel.h"
#include "polynomialBasis.h"
#include "JacobianBasis.h"
#include "AnalyseCurvedMesh.h"
#include "fullMatrix.h"
extern "C"
{
......@@ -24,26 +26,260 @@ extern "C"
std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
{
return "Plugin(AnalyseCurvedMesh) computes AnalyseCurvedMeshs to boundaries in "
"a mesh.\n\n"
"Plugin(AnalyseCurvedMesh) creates a bunch of files.";
return "Plugin(AnalyseCurvedMesh) check the jacobian of all elements of the greater dimension. Elements for which we are absolutely certain that they are good (positive jacobian) are ignored. Others are only suppose to be wrong."
"Plugin(AnalyseCurvedMesh) write in the console which elements are supposed to be wrong. (if labels of analysed type of elements are set visible, only supposed wrong elements will be visible)";
}
static double computeDeterminant(MElement *el, double jac[3][3])
{
switch (el->getDim()) {
case 0:
return 1.0;
case 1:
return jac[0][0];
case 2:
return jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
case 3:
return jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2];
default:
return 1;
}
}
double getJacobian(double gsf[][3], double jac[3][3], MElement *el)
{
jac[0][0] = jac[0][1] = jac[0][2] = 0.;
jac[1][0] = jac[1][1] = jac[1][2] = 0.;
jac[2][0] = jac[2][1] = jac[2][2] = 0.;
for (int i = 0; i < el->getNumVertices(); i++) {
const MVertex *v = el->getVertex(i);
double *gg = gsf[i];
for (int j = 0; j < 3; j++) {
jac[j][0] += v->x() * gg[j];
jac[j][1] += v->y() * gg[j];
jac[j][2] += v->z() * gg[j];
}
}
return computeDeterminant(el, jac);
}
PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
{
Msg::Info("Hello !");
Msg::Info("AnalyseCurvedMeshPlugin : Starting analyse.");
int numBadEl = 0;
int numAnalysedEl = 0;
GModel *m = GModel::current();
JacobianBases::find(MSH_QUA_8);
JacobianBases::find(MSH_QUA_9);
JacobianBases::find(MSH_TET_20);
JacobianBases::find(MSH_PRI_18);
switch (m->getDim()) {
case 3 :
{
Msg::Info("Only 3D elements will be analyse.");
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++) {
GRegion *r = *it;
unsigned int numType[5] = {0, 0, 0, 0, 0};
r->getNumMeshElements(numType);
for(int type = 0; type < 5; type++) { // loop over all types of elements
MElement *const *el = r->getStartElementType(type);
for(unsigned int i = 0; i < numType[type]; i++) { // loop over all elements of a type
numAnalysedEl++;
if (checkJacobian((*(el+i)), 0) <= 0) numBadEl++;
}
}
}
break;
}
case 2 :
{
Msg::Info("Only 2D elements will be analyse.");
Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) {
GFace *f = *it;
Msg::Info("Je suis dans la surface %d", f->tag());
for(int i = 0; i < (int) f->triangles.size(); i++) {
MTriangle *t = f->triangles[i];
const polynomialBasis *p = t->getFunctionSpace();
Msg::Info("Taille de la matrice P : %dx%d", p->monomials.size1(), p->monomials.size2());
unsigned int numType[3] = {0, 0, 0};
f->getNumMeshElements(numType);
for (int type = 0; type < 3; type++) {
MElement *const *el = f->getStartElementType(type);
for (unsigned int i = 0; i < numType[type]; i++) {
numAnalysedEl++;
if (checkJacobian((*(el+i)), 0) <= 0) numBadEl++;
}
}
}
break;
}
case 1 :
{
Msg::Info("Only 1D elements will be analyse.");
Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) {
GEdge *e = *it;
unsigned int *numElement = {0};
e->getNumMeshElements(numElement);
MElement *const *el = e->getStartElementType(0);
for (unsigned int i = 0; i < *numElement; i++) {
numAnalysedEl++;
if (checkJacobian((*(el+i)), 0) <= 0) numBadEl++;
}
}
break;
}
default :
{
Msg::Error("I can't analyse any element.");
}
}
Msg::Info("%d elements have been analysed.", numAnalysedEl);
Msg::Info("%d elements may be bad.", numBadEl);
Msg::Info("AnalyseCurvedMeshPlugin : Job finished.");
return 0;
}
bool GMSH_AnalyseCurvedMeshPlugin::isJacPositive(MElement *el)
{
const polynomialBasis *fs = el->getFunctionSpace(-1);
if (!fs) {
Msg::Error("Function space not implemented for type of element %d", el->getNum());
return false;
}
const JacobianBasis *jfs = el->getJacobianFuncSpace(-1);
if (!jfs) {
Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum());
return false;
}
int numSamplingPt = jfs->points.size1();
int dim = jfs->points.size2();
bool isPositive = true;
fullVector<double> jacobian(numSamplingPt);
for (int i = 0; i < numSamplingPt; i++) {
double gsf[256][3];
switch (dim) {
case 1:
fs->df(jfs->points(i,0),0,0, gsf);
break;
case 2:
fs->df(jfs->points(i,0),jfs->points(i,1),0, gsf);
break;
case 3:
fs->df(jfs->points(i,0),jfs->points(i,1),jfs->points(i,2), gsf);
break;
default:
Msg::Error("Can't get the gradient for %dD elements.", dim);
return false;
}
double jac[3][3];
jacobian(i) = getJacobian(gsf, jac, el);
}
fullVector<double> jacBez(jacobian.size());
jfs->matrixLag2Bez.mult(jacobian, jacBez);
for (int i = 0; i < jacBez.size(); i++) {
if (jacBez(i) < 0.) isPositive = false;
}
//Msg::Info("%d", el->getNum());
//if (el->getNum() == 582 || el->getNum() == 584) {
// for (int i = 0; i < jfs->points.size1(); i++) {
// Msg::Info("jacBez[%d] = %lf, jacobian[%d] = %lf", i, jacBez(i), i, jacobian(i));
// }
//}
return isPositive;
}
int GMSH_AnalyseCurvedMeshPlugin::checkJacobian(MElement *el, int depth)
{
int tag = method1(el, depth);
if (tag < 0) {
Msg::Info("Bad element : %d", el->getNum());
}
else if (tag > 0) {
el->setVisibility(0);
}
else {
Msg::Info("Element %d may be bad.", el->getNum());
}
return tag;
}
int GMSH_AnalyseCurvedMeshPlugin::method1(MElement *el, int depth)
{
const polynomialBasis *fs = el->getFunctionSpace(-1);
if (!fs) {
Msg::Error("Function space not implemented for type of element %d", el->getNum());
return 0;
}
const JacobianBasis *jfs = el->getJacobianFuncSpace(-1);
if (!jfs) {
Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum());
return 0;
}
int numSamplingPt = jfs->points.size1(), dim = jfs->points.size2();
fullVector<double> jacobian(numSamplingPt);
for (int i = 0; i < numSamplingPt; i++) {
double gsf[256][3];
switch (dim) {
case 1:
fs->df(jfs->points(i,0),0,0, gsf);
break;
case 2:
fs->df(jfs->points(i,0),jfs->points(i,1),0, gsf);
break;
case 3:
fs->df(jfs->points(i,0),jfs->points(i,1),jfs->points(i,2), gsf);
break;
default:
Msg::Error("Can't get the gradient for %dD elements.", dim);
return false;
}
double jac[3][3];
jacobian(i) = getJacobian(gsf, jac, el);
}
for (int i = 0; i < jacobian.size(); i++) {
if (jacobian(i) < 0.) return -1;
}
fullVector<double> jacBez(jacobian.size());
jfs->matrixLag2Bez.mult(jacobian, jacBez);
bool allPtPositive = true;
for (int i = 0; i < jacBez.size(); i++) {
if (jacBez(i) < 0.) allPtPositive = false;
}
if (allPtPositive) return 1;
if (depth <= 0) {
return 0;
}
//else{return division(el, depth-1);}
return 0;
}
\ No newline at end of file
Plugin/AnalyseCurvedMesh.h
+ 10
− 1
View file @ 4d2b0ea1
......@@ -3,13 +3,19 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
//
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
// Contributed by Amaury J.
// Added code :
// Numeric/polynomialBasis.cpp -> jacobianPolynomialBases::find(int tag), #include
// Geo/MElement.h -> getJacobianFunctionSpace(int), #include
// Geo/MTriangle .h/.cpp -> getJacobianFunctionSpace(int)
#ifndef _ANALYSECURVEDMESH_H_
#define _ANALYSECURVEDMESH_H_
#include <string>
#include "Plugin.h"
#include "MElement.h"
extern "C"
{
......@@ -28,6 +34,9 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
std::string getHelp() const;
std::string getAuthor() const { return "Amaury Johnen"; }
PView *execute(PView *);
bool isJacPositive(MElement *);
int method1(MElement *, int depth);
int checkJacobian(MElement *, int depth);
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment