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

AnalyseCurvedMesh development (should work in all cases, edit...

AnalyseCurvedMesh development (should work in all cases, edit PluginManager.cpp to activate the plugin)
parent d3c40288
No related branches found
No related tags found
No related merge requests found
......@@ -119,6 +119,11 @@ class MPolyhedron : public MElement {
if(_orig) return _orig->getFunctionSpace(order);
return 0;
}
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
{
if(_orig) return _orig->getJacobianFuncSpace(order);
return 0;
}
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
......@@ -241,6 +246,11 @@ class MPolygon : public MElement {
if(_orig) return _orig->getFunctionSpace(order);
return 0;
}
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
{
if(_orig) return _orig->getJacobianFuncSpace(order);
return 0;
}
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
......@@ -297,6 +307,11 @@ class MLineChild : public MLine {
if(_orig) return _orig->getFunctionSpace(order);
return 0;
}
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
{
if(_orig) return _orig->getJacobianFuncSpace(order);
return 0;
}
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
......
......@@ -29,6 +29,26 @@ const polynomialBasis* MLine::getFunctionSpace(int o) const
return 0;
}
const JacobianBasis* MLine::getJacobianFuncSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
case 1: return JacobianBases::find(MSH_LIN_2);
case 2: return JacobianBases::find(MSH_LIN_3);
case 3: return JacobianBases::find(MSH_LIN_4);
case 4: return JacobianBases::find(MSH_LIN_5);
case 5: return JacobianBases::find(MSH_LIN_6);
case 6: return JacobianBases::find(MSH_LIN_7);
case 7: return JacobianBases::find(MSH_LIN_8);
case 8: return JacobianBases::find(MSH_LIN_9);
case 9: return JacobianBases::find(MSH_LIN_10);
case 10: return JacobianBases::find(MSH_LIN_11);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
}
void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = getNGQLPts(pOrder);
......
......@@ -70,6 +70,7 @@ class MLine : public MElement {
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = 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;
......
......@@ -52,6 +52,7 @@ class MPoint : public MElement {
}
virtual const polynomialBasis* getFunctionSpace(int o) const { return polynomialBases::find(MSH_PNT); }
virtual const JacobianBasis* getJacobianFuncSpace(int o) const { return JacobianBases::find(MSH_PNT); }
virtual bool isInside(double u, double v, double w)
{
return true;
......
......@@ -53,6 +53,29 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
return 0;
}
const JacobianBasis* MPrism::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_PRI_6);
case 2: return JacobianBases::find(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
else {
switch (order) {
case 1: return JacobianBases::find(MSH_PRI_6);
case 2: return JacobianBases::find(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
return 0;
}
double MPrism::getInnerRadius()
{
double dist[3], k = 0.;
......
......@@ -133,6 +133,7 @@ class MPrism : public MElement {
tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual int getVolumeSign();
/* virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
......
......@@ -52,6 +52,42 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
return 0;
}
const JacobianBasis* MQuadrangle::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_QUA_4);
case 2: return JacobianBases::find(MSH_QUA_8);
case 3: return JacobianBases::find(MSH_QUA_12);
case 4: return JacobianBases::find(MSH_QUA_16I);
case 5: return JacobianBases::find(MSH_QUA_20);
case 6: return JacobianBases::find(MSH_QUA_24);
case 7: return JacobianBases::find(MSH_QUA_28);
case 8: return JacobianBases::find(MSH_QUA_32);
case 9: return JacobianBases::find(MSH_QUA_36I);
case 10: return JacobianBases::find(MSH_QUA_40);
}
}
switch (order) {
case 1: return JacobianBases::find(MSH_QUA_4);
case 2: return JacobianBases::find(MSH_QUA_9);
case 3: return JacobianBases::find(MSH_QUA_16);
case 4: return JacobianBases::find(MSH_QUA_25);
case 5: return JacobianBases::find(MSH_QUA_36);
case 6: return JacobianBases::find(MSH_QUA_49);
case 7: return JacobianBases::find(MSH_QUA_64);
case 8: return JacobianBases::find(MSH_QUA_81);
case 9: return JacobianBases::find(MSH_QUA_100);
case 10: return JacobianBases::find(MSH_QUA_121);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
return 0;
}
int MQuadrangleN::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
int MQuadrangle8::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
int MQuadrangle9::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
......
......@@ -119,6 +119,7 @@ class MQuadrangle : public MElement {
virtual const char *getStringForDIFF() const { return "ElmB4n2D"; }
virtual const char *getStringForINP() const { return "C2D4"; }
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
......
This diff is collapsed.
......@@ -6,21 +6,22 @@
#ifndef _JACOBIAN_BASIS_H_
#define _JACOBIAN_BASIS_H_
#include <math.h>
#include <map>
#include <vector>
#include "fullMatrix.h"
#include <iostream>
class JacobianBasis
{
public:
int numLagPts;
int numDivisions;
fullMatrix<double> exposants; //exposants of Bezier FF
fullMatrix<double> points; //sampling point
fullMatrix<double> matrixLag2Bez;
fullMatrix<double> gradShapeMatX;
fullMatrix<double> gradShapeMatY;
fullMatrix<double> gradShapeMatZ;
fullMatrix<double> divisor;
};
class JacobianBases
......@@ -31,207 +32,4 @@ class JacobianBases
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
......@@ -42,7 +42,7 @@ extern "C" {
}
template<>
void fullVector<double>::axpy(fullVector<double> &x,double alpha)
void fullVector<double>::axpy(const fullVector<double> &x,double alpha)
{
int M = _r, INCX = 1, INCY = 1;
F77NAME(daxpy)(&M, &alpha, x._data,&INCX, _data, &INCY);
......@@ -220,7 +220,7 @@ bool fullMatrix<double>::luSolve(const fullVector<double> &rhs, fullVector<doubl
}
template<>
bool fullMatrix<double>::invert(fullMatrix<double> &result)
bool fullMatrix<double>::invert(fullMatrix<double> &result) const
{
int M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)];
......
......@@ -47,6 +47,13 @@ class fullVector
}
return false;
}
void setAsProxy(const fullVector<scalar> &original, int r_start, int r)
{
if(_data)
delete [] _data;
_r = r;
_data = original._data + r_start;
}
inline scalar operator () (int i) const
{
return _data[i];
......@@ -77,13 +84,17 @@ class fullVector
for(int i = 0; i < _r; ++i) s += _data[i] * other._data[i];
return s;
}
void axpy(fullVector<scalar> &x, scalar alpha=1.)
void axpy(const fullVector<scalar> &x, scalar alpha=1.)
#if !defined(HAVE_BLAS)
{
for (int i = 0; i < _r; i++) _data[i] += alpha * x._data[i];
}
#endif
;
void multTByT(const fullVector<scalar> &x)
{
for (int i = 0; i < _r; i++) _data[i] *= x._data[i];
}
void print(const char *name="") const
{
printf("Printing vector %s:\n", name);
......@@ -260,6 +271,16 @@ class fullMatrix
for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
(*this)(desti, destj) = a(i, j);
}
void copy(const fullMatrix<scalar> &a)
{
if(_data && _own_data)
delete [] _data;
_r = a._r;
_c = a._c;
_data = new scalar[_r * _c];
_own_data = true;
setAll(a);
}
void mult_naive(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)const
{
c.scale(0.);
......@@ -276,6 +297,10 @@ class fullMatrix
}
#endif
;
void multTByT(const fullMatrix<scalar> &a)
{
for (int i = 0; i < _r * _c; i++) _data[i] *= a._data[i];
}
void gemm_naive(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
scalar alpha=1., scalar beta=1.)
{
......@@ -386,7 +411,7 @@ class fullMatrix
}
#endif
;
bool invert(fullMatrix<scalar> &result)
bool invert(fullMatrix<scalar> &result) const
#if !defined(HAVE_LAPACK)
{
Msg::Error("LU factorization requires LAPACK");
......
This diff is collapsed.
......@@ -2,20 +2,13 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
//
// 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"
#include <vector>
extern "C"
{
......@@ -33,10 +26,19 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
}
std::string getHelp() const;
std::string getAuthor() const { return "Amaury Johnen"; }
int getNbOptions() const;
StringXNumber *getOption(int);
PView *execute(PView *);
bool isJacPositive(MElement *);
int method1(MElement *, int depth);
int checkJacobian(MElement *, int depth);
int method_1_1(MElement *, int depth);
int method_1_2(MElement *, int depth);
int method_1_3(MElement *, int depth);
void method_2_2(MElement *const *, std::vector<int> &tags, int depth);
void method_2_3(MElement *const *, std::vector<int> &tags, int depth);
//int checkJacobian(MElement *, int depth);
//int *checkJacobian2(MElement *const *, int numEl, int depth);
int *checkJacobian(MElement *const *, int numEl, int depth, int method);
int division(const JacobianBasis *, const fullVector<double> &, int depth);
};
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment