Commit 201ba8a5 authored by Christophe Geuzaine's avatar Christophe Geuzaine

Merge branch 'curveBL' into 'master'

Update of boundary layer curver

See merge request !141
parents 190b4ff3 f2920fed
Pipeline #1837 passed with stage
in 62 minutes and 7 seconds
......@@ -9,4 +9,8 @@ contrib/Parasolid
*~
*.so
*.so.*
*#
\ No newline at end of file
*#
## CLion directories
.idea
cmake-build-*
......@@ -221,8 +221,12 @@ static void highordertools_runopti_cb(Fl_Widget *w, void *data)
case 3: { // Fast curving 2
FastCurvingParameters p;
p.onlyVisible = onlyVisible;
p.dim = dim;
p.thickness = true;
if (dim == 3) {
p.dim = 2;
HighOrderMeshFastCurving(GModel::current(), p);
}
p.dim = dim;
HighOrderMeshFastCurving(GModel::current(), p);
break;
}
......@@ -333,7 +337,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
{"Optimization", 0, 0, 0},
{"Elastic Analogy", 0, 0, 0},
{"Fast Curving", 0, 0, 0},
{"Fast Curving 2 (experimental)", 0, 0, 0},
{"Curve boundary layers (3D is experimental)", 0, 0, 0},
{0}
};
choice[2] = new Fl_Choice
......
......@@ -6,6 +6,10 @@
#include <algorithm>
#include "MEdge.h"
#include "Numeric.h"
#include "GmshDefines.h"
#include "ElementType.h"
#include "nodalBasis.h"
#include "BasisFactory.h"
// FIXME
// remove that when MElementCut is removed
......@@ -64,6 +68,7 @@ bool MEdge::isInside(MVertex *v) const
return true;
}
bool SortEdgeConsecutive(const std::vector<MEdge> &e,
std::vector<std::vector<MVertex*> >&vs)
{
......@@ -87,10 +92,10 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
std::map<MVertex*, std::pair<MVertex*,MVertex*> >::iterator it = c.begin();
start = it->first;
for (; it != c.end(); ++it) {
if (it->second.second == NULL){
start = it->first;
break;
}
if (it->second.second == NULL){
start = it->first;
break;
}
}
}
std::map<MVertex*, std::pair<MVertex*,MVertex*> >::iterator it = c.find(start);
......@@ -106,11 +111,11 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
if (v1 == prev)current = v2;
else if (v2 == prev)current = v1;
else {
break;
break;
}
prev = temp;
if (current == start) {
v.push_back(current);
v.push_back(current);
}
} while (current != start && current != NULL);
vs.push_back(v);
......@@ -118,3 +123,64 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
return true;
}
MEdgeN::MEdgeN(const std::vector<MVertex*> &v)
{
_v.resize(v.size());
for(unsigned int i = 0; i < v.size(); i++)
_v[i] = v[i];
}
MEdge MEdgeN::getEdge() const
{
return MEdge(_v[0], _v[1]);
}
SPoint3 MEdgeN::pnt(double u) const
{
int tagLine = ElementType::getType(TYPE_LIN, getPolynomialOrder());
const nodalBasis *fs = BasisFactory::getNodalBasis(tagLine);
double f[100];
fs->f(u, 0, 0, f);
double x = 0, y = 0, z = 0;
for (int i = 0; i < fs->getNumShapeFunctions(); i++) {
x += f[i] * _v[i]->x();
y += f[i] * _v[i]->y();
z += f[i] * _v[i]->z();
}
return SPoint3(x, y, z);
}
SVector3 MEdgeN::tangent(double u) const
{
int tagLine = ElementType::getType(TYPE_LIN, getPolynomialOrder());
const nodalBasis *fs = BasisFactory::getNodalBasis(tagLine);
double sf[100][3];
fs->df(u, 0, 0, sf);
double dx = 0, dy = 0, dz = 0;
for (int i = 0; i < fs->getNumShapeFunctions(); i++) {
dx += sf[i][0] * _v[i]->x();
dy += sf[i][0] * _v[i]->y();
dz += sf[i][0] * _v[i]->z();
}
return SVector3(dx, dy, dz).unit();
}
double MEdgeN::interpolate(const double val[], double u, int stride) const
{
int tagLine = ElementType::getType(TYPE_LIN, getPolynomialOrder());
const nodalBasis *fs = BasisFactory::getNodalBasis(tagLine);
double f[100];
fs->f(u, 0, 0, f);
double sum = 0;
for (int i = 0, k = 0; i < fs->getNumShapeFunctions(); i++, k += stride) {
sum += f[i] * val[k];
}
return sum;
}
......@@ -127,4 +127,24 @@ struct Less_Edge : public std::binary_function<MEdge, MEdge, bool> {
bool SortEdgeConsecutive(const std::vector<MEdge> &,
std::vector<std::vector<MVertex*> >&vs);
class MEdgeN {
private:
std::vector<MVertex *> _v;
public:
MEdgeN() {}
MEdgeN(const std::vector<MVertex*> &v);
inline int getNumVertices() const { return (int)_v.size(); }
inline MVertex *getVertex(int i) const { return _v[i]; }
inline const std::vector<MVertex*> &getVertices() const { return _v; }
inline int getPolynomialOrder() const {return getNumVertices() - 1;}
MEdge getEdge() const;
SPoint3 pnt(double u) const;
SVector3 tangent(double u) const;
double interpolate(const double val[], double u, int stride = 1) const;
};
#endif
......@@ -75,6 +75,46 @@ double MElement::getTolerance()
return _isInsideTolerance;
}
bool MElement::_getFaceInfo(const MFace &face, const MFace &other,
int &sign, int &rot)
{
// Looks how is 'other' compared to 'face'. We suppose that 'face'
// is the reference.
// Return false if they are different.
// In case sign = 1, then rot = 1 if 'other' is equal to 'face' rotated
// one time according to the right hand side.
// In case sign = -1, then rot = 0 if 'other' is equal to reversed 'face'.
// In case sign = -1, then rot = 1 if 'other' is equal to reversed 'face'
// rotated one time according to the right hand side.
int N = face.getNumVertices();
sign = 0;
rot = -1;
if (N != other.getNumVertices()) return false;
sign = 1;
for (rot = 0; rot < N; ++rot) {
int i;
for (i = 0; i < N; ++i) {
if (other.getVertex(i) != face.getVertex((i+rot)%N)) break;
}
if (i == N) return true;
}
sign = -1;
for (rot = 0; rot < N; ++rot) {
int i;
for (i = 0; i < N; ++i) {
if (other.getVertex(i) != face.getVertex((N+rot-i)%N)) break;
}
if (i == N) return true;
}
sign = 0;
rot = -1;
return false;
}
void MElement::_getEdgeRep(MVertex *v0, MVertex *v1,
double *x, double *y, double *z, SVector3 *n,
int faceIndex)
......@@ -124,6 +164,73 @@ char MElement::getVisibility() const
return _visible;
}
MEdgeN MElement::getHighOrderEdge(int num, int sign)
{
const int order = getPolynomialOrder();
std::vector<MVertex*> vertices((unsigned int) order + 1);
vertices[0] = getVertex(numEdge2numVertex(num, sign > 0 ? 0 : 1));
vertices[1] = getVertex(numEdge2numVertex(num, sign > 0 ? 1 : 0));
const int start = getNumPrimaryVertices() + num * (order - 1);
const int end = getNumPrimaryVertices() + (num + 1) * (order - 1);
int k = 1;
if (sign > 0) {
for (int i = start; i < end; ++i) {
vertices[++k] = getVertex(i);
}
}
else {
for (int i = end - 1; i >= start; --i) {
vertices[++k] = getVertex(i);
}
}
return MEdgeN(vertices);
}
bool MElement::getEdgeInfo(const MEdge &edge, int &ithEdge, int &sign) const
{
for (ithEdge = 0; ithEdge < getNumEdges(); ithEdge++) {
const MVertex *v0 = getVertex(numEdge2numVertex(ithEdge, 0));
const MVertex *v1 = getVertex(numEdge2numVertex(ithEdge, 1));
if (v0 == edge.getVertex(0) && v1 == edge.getVertex(1)){
sign = 1; return true;
}
if (v1 == edge.getVertex(0) && v0 == edge.getVertex(1)){
sign = -1; return true;
}
}
Msg::Error("Could not get edge information for element %d", getNum());
return false;
}
MFaceN MElement::getHighOrderFace(int num, int sign, int rot)
{
if (getDim() < 2 || getDim() > 3) {
Msg::Error("Wrong dimension for getHighOrderFace");
return MFaceN();
}
if (getDim() == 2) {
std::vector<MVertex *> vertices(getNumVertices());
getVertices(vertices);
return MFaceN(getType(), getPolynomialOrder(), vertices);
}
const nodalBasis *fs = getFunctionSpace();
int id = fs->getClosureId(num, sign, rot);
const std::vector<int> &closure = fs->getClosure(id);
std::vector<MVertex *> vertices(closure.size());
for (unsigned int i = 0; i < closure.size(); ++i) {
vertices[i] = getVertex(closure[i]);
}
static int type2numTriFaces[9] = {0, 0, 0, 1, 0, 4, 4, 2, 0};
int typeFace = TYPE_TRI;
if (num >= type2numTriFaces[getType()]) typeFace = TYPE_QUA;
return MFaceN(typeFace, getPolynomialOrder(), vertices);
}
double MElement::minEdge()
{
double m = 1.e25;
......@@ -362,15 +469,12 @@ void MElement::signedInvGradErrorRange(double &minSIGE, double &maxSIGE)
void MElement::getNode(int num, double &u, double &v, double &w) const
{
// only for MElements that don't have a lookup table for this
// (currently only 1st order elements have)
double uvw[3];
const MVertex* ver = getVertex(num);
double xyz[3] = {ver->x(), ver->y(), ver->z()};
xyz2uvw(xyz, uvw);
u = uvw[0];
v = uvw[1];
w = uvw[2];
// Should we always do this instead of using lookup table for linear elements?
const nodalBasis *nb = getFunctionSpace();
const fullMatrix<double> &refpnts = nb->getReferenceNodes();
u = refpnts(num, 0);
v = getDim() > 1 ? refpnts(num, 1) : 0;
v = getDim() > 2 ? refpnts(num, 2) : 0;
}
void MElement::getShapeFunctions(double u, double v, double w, double s[], int o) const
......@@ -578,7 +682,8 @@ const nodalBasis* MElement::getFunctionSpace(int order, bool serendip) const
const JacobianBasis* MElement::getJacobianFuncSpace(int order) const
{
if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
int tag = ElementType::getType(getType(), order);
return BasisFactory::getJacobianBasis(tag);
}
static double _computeDeterminantAndRegularize(const MElement *ele, double jac[3][3])
......@@ -999,12 +1104,15 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3]) const
}
double inv[3][3];
inv3x3(jac, inv);
double un = uvw[0] + inv[0][0] * (xyz[0] - xn) +
inv[1][0] * (xyz[1] - yn) + inv[2][0] * (xyz[2] - zn);
double vn = uvw[1] + inv[0][1] * (xyz[0] - xn) +
inv[1][1] * (xyz[1] - yn) + inv[2][1] * (xyz[2] - zn);
double wn = uvw[2] + inv[0][2] * (xyz[0] - xn) +
inv[1][2] * (xyz[1] - yn) + inv[2][2] * (xyz[2] - zn);
double un = uvw[0] + inv[0][0] * (xyz[0] - xn)
+ inv[1][0] * (xyz[1] - yn)
+ inv[2][0] * (xyz[2] - zn);
double vn = uvw[1] + inv[0][1] * (xyz[0] - xn)
+ inv[1][1] * (xyz[1] - yn)
+ inv[2][1] * (xyz[2] - zn);
double wn = uvw[2] + inv[0][2] * (xyz[0] - xn)
+ inv[1][2] * (xyz[1] - yn)
+ inv[2][2] * (xyz[2] - zn);
error = sqrt(SQU(un - uvw[0]) + SQU(vn - uvw[1]) + SQU(wn - uvw[2]));
uvw[0] = un;
uvw[1] = vn;
......
......@@ -49,6 +49,10 @@ class MElement
void _getFaceRepQuad(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
double *x, double *y, double *z, SVector3 *n);
#endif
static bool _getFaceInfo(const MFace &face, const MFace &other,
int &sign, int &rot);
public :
MElement(int num=0, int part=0);
virtual ~MElement(){}
......@@ -147,11 +151,20 @@ class MElement
// get the edges
virtual int getNumEdges() const = 0;
virtual MEdge getEdge(int num) const= 0;
virtual MEdgeN getHighOrderEdge(int num, int sign);
virtual MEdgeN getHighOrderEdge(const MEdge &edge)
{
int num, sign;
if (!getEdgeInfo(edge, num, sign)) return MEdgeN();
return getHighOrderEdge(num, sign);
}
// give an MEdge as input and get its local number and sign
virtual void getEdgeInfo(const MEdge & edge, int &ithEdge, int &sign) const
virtual bool getEdgeInfo(const MEdge & edge, int &ithEdge, int &sign) const;
virtual int numEdge2numVertex(int numEdge, int numVert) const
{
Msg::Error("Edge information not available for this element");
return -1;
}
// get an edge representation for drawing
......@@ -161,13 +174,22 @@ class MElement
// get the faces
virtual int getNumFaces() = 0;
virtual MFace getFace(int num) = 0;
virtual MFace getFace(int num) const = 0;
virtual MFaceN getHighOrderFace(int num, int sign, int rot);
virtual MFaceN getHighOrderFace(const MFace &face)
{
int num, sign, rot;
if (!getFaceInfo(face, num, sign, rot))
return MFaceN();
return this->getHighOrderFace(num, sign, rot);
}
// give an MFace as input and get its local number, sign and rotation
virtual void getFaceInfo(const MFace & face, int &ithFace, int &sign,
virtual bool getFaceInfo(const MFace & face, int &ithFace, int &sign,
int &rot) const
{
Msg::Error("Face information not available for this element");
return false;
}
// get a face representation for drawing
......
......@@ -88,7 +88,7 @@ class MPolyhedron : public MElement {
v[1] = _edges[num].getVertex(1);
}
virtual int getNumFaces() { return _faces.size(); }
virtual MFace getFace(int num) { return _faces[num]; }
virtual MFace getFace(int num) const { return _faces[num]; }
virtual int getNumFacesRep(bool curved) { return _faces.size(); }
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
{
......@@ -246,7 +246,7 @@ class MPolygon : public MElement {
v[1] = _edges[num].getVertex(1);
}
virtual int getNumFaces() { return 1; }
virtual MFace getFace(int num) { return MFace(_vertices); }
virtual MFace getFace(int num) const { return MFace(_vertices); }
virtual int getNumFacesRep(bool curved) { return _parts.size(); }
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
{
......
......@@ -8,6 +8,9 @@
#include "GmshConfig.h"
#include "MFace.h"
#include "Numeric.h"
#include "ElementType.h"
#include "nodalBasis.h"
#include "BasisFactory.h"
bool compare(const MVertex *const v0, const MVertex *const v1){ return v0->getNum() < v1->getNum(); }
......@@ -50,7 +53,7 @@ MFace::MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3)
sortVertices(_v, _si);
}
MFace::MFace(std::vector<MVertex*> v)
MFace::MFace(const std::vector<MVertex*> &v)
{
for(unsigned int i = 0; i < v.size(); i++)
_v.push_back(v[i]);
......@@ -95,3 +98,173 @@ bool MFace::computeCorrespondence(const MFace &other,
}
return false;
}
MFaceN::MFaceN(int type, int order, const std::vector<MVertex*> &v)
: _type(type), _order(order)
{
_v.resize(v.size());
for(unsigned int i = 0; i < v.size(); i++)
_v[i] = v[i];
}
MEdgeN MFaceN::getHighOrderEdge(int num, int sign) const
{
int nCorner = getNumCorners();
std::vector<MVertex*> vertices((unsigned int)_order + 1);
if (sign == 1) {
vertices[0] = _v[num];
vertices[1] = _v[(num + 1) % nCorner];
}
else {
vertices[0] = _v[(num + 1) % nCorner];
vertices[1] = _v[num];
}
int start = nCorner + num * (_order - 1);
int end = nCorner + (num + 1) * (_order - 1);
int k = 1;
if (sign == 1) {
for (int i = start; i < end; ++i) vertices[++k] = _v[i];
}
else {
for (int i = end-1; i >= start; --i) vertices[++k] = _v[i];
}
return MEdgeN(vertices);
}
MFace MFaceN::getFace() const
{
if (_type == TYPE_TRI)
return MFace(_v[0], _v[1], _v[2]);
else
return MFace(_v[0], _v[1], _v[2], _v[3]);
}
SPoint3 MFaceN::pnt(double u, double v) const
{
int tag = ElementType::getType(_type, _order);
const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
double f[100];
fs->f(u, v, 0, f);
double x = 0, y = 0, z = 0;
for (int j = 0; j < fs->getNumShapeFunctions(); j++) {
x += f[j] * _v[j]->x();
y += f[j] * _v[j]->y();
z += f[j] * _v[j]->z();
}
return SPoint3(x, y, z);
}
SVector3 MFaceN::tangent(double u, double v, int num) const
{
if (num != 0 && num != 1) num = 0;
int tag = ElementType::getType(_type, _order);
const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
double sf[100][3];
fs->df(u, v, 0, sf);
double dx = 0, dy = 0, dz = 0;
for (int j = 0; j < fs->getNumShapeFunctions(); j++) {
dx += sf[j][num] * _v[j]->x();
dy += sf[j][num] * _v[j]->y();
dz += sf[j][num] * _v[j]->z();
}
return SVector3(dx, dy, dz).unit();
}
SVector3 MFaceN::normal(double u, double v) const
{
int tag = ElementType::getType(_type, _order);
const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
double sf[100][3];
fs->df(u, v, 0, sf);
double dx[2] = {0, 0}, dy[2] = {0, 0}, dz[2] = {0, 0};
for (int j = 0; j < fs->getNumShapeFunctions(); j++) {
for (int k = 0; k < 1; ++k) {
dx[k] += sf[j][k] * _v[j]->x();
dy[k] += sf[j][k] * _v[j]->y();
dz[k] += sf[j][k] * _v[j]->z();
}
}
SVector3 t0 = SVector3(dx[0], dy[0], dz[0]);
SVector3 t1 = SVector3(dx[1], dy[1], dz[1]);
return crossprod(t0, t1).unit();
}
void MFaceN::frame(double u, double v,
SVector3 &t0, SVector3 &t1, SVector3 &n) const
{
int tag = ElementType::getType(_type, _order);
const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
double sf[100][3];
fs->df(u, v, 0, sf);
double dx[2] = {0, 0}, dy[2] = {0, 0}, dz[2] = {0, 0};
for (int j = 0; j < fs->getNumShapeFunctions(); j++) {
for (int k = 0; k < 2; ++k) {
dx[k] += sf[j][k] * _v[j]->x();
dy[k] += sf[j][k] * _v[j]->y();
dz[k] += sf[j][k] * _v[j]->z();
}
}
t0 = SVector3(dx[0], dy[0], dz[0]).unit();
t1 = SVector3(dx[1], dy[1], dz[1]).unit();
n = crossprod(t0, t1);
}
void MFaceN::frame(double u, double v, SPoint3 &p,
SVector3 &t0, SVector3 &t1, SVector3 &n) const
{
int tag = ElementType::getType(_type, _order);
const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
double f[100];
double sf[100][3];
fs->f(u, v, 0, f);
fs->df(u, v, 0, sf);
double x = 0, y = 0, z = 0;
double dx[2] = {0, 0}, dy[2] = {0, 0}, dz[2] = {0, 0};
for (int j = 0; j < fs->getNumShapeFunctions(); j++) {
x += f[j] * _v[j]->x();
y += f[j] * _v[j]->y();
z += f[j] * _v[j]->z();
for (int k = 0; k < 2; ++k) {
dx[k] += sf[j][k] * _v[j]->x();
dy[k] += sf[j][k] * _v[j]->y();
dz[k] += sf[j][k] * _v[j]->z();
}
}
p = SPoint3(x, y, z);
t0 = SVector3(dx[0], dy[0], dz[0]).unit();
t1 = SVector3(dx[1], dy[1], dz[1]).unit();
n = crossprod(t0, t1);
}
void MFaceN::repositionInnerVertices(const fullMatrix<double> *placement) const
{
int nCorner = getNumCorners();
int start = nCorner + (_order - 1) * nCorner;
for (int i = start; i < (int)_v.size(); ++i) {
MVertex *v = _v[i];
v->x() = 0;
v->y() = 0;
v->z() = 0;
for (int j = 0; j < placement->size2(); ++j) {
const double coeff = (*placement)(i - start, j);
v->x() += coeff * _v[j]->x();
v->y() += coeff * _v[j]->y();
v->z() += coeff * _v[j]->z();
}
}
}
......@@ -12,6 +12,10 @@
#include "MEdge.h"
#include "SVector3.h"
#include "GmshMessage.h"
#include "GmshDefines.h"
template <class t> class fullMatrix;
// A mesh face.
class MFace {
......@@ -22,7 +26,7 @@ class MFace {
public:
MFace() {}
MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3=0);
MFace(const std::vector<MVertex*> v);
MFace(const std::vector<MVertex*> &v);
inline int getNumVertices() const { return (int)_v.size(); }
inline MVertex *getVertex(const int i) const { return _v[i]; }
inline MVertex *getSortedVertex(const int i) const { return _v[int(_si[i])]; }
......@@ -142,5 +146,37 @@ struct Less_Face : public std::binary_function<MFace, MFace, bool> {
}
};
class MFaceN {
private:
int _type;
int _order;
std::vector<MVertex *> _v;
public:
MFaceN() {}
MFaceN(int type, int order, const std::vector<MVertex*> &v);
inline int getPolynomialOrder() const {return _order;}
inline int getType() const {return _type;}
inline bool isTriangular() const {return _type == TYPE_TRI;}
inline int getNumVertices() const { return (int)_v.size(); }
inline int getNumCorners() const { return isTriangular() ? 3 : 4; }
inline int getNumVerticesOnBoundary() const { return getNumCorners() * _order; }
inline MVertex *getVertex(int i) const { return _v[i]; }
inline const std::vector<MVertex*> &getVertices() const { return _v; }
MEdgeN getHighOrderEdge(int num, int sign) const;
MFace getFace() const;
SPoint3 pnt(double u, double v) const;
SVector3 tangent(double u, double v, int num) const;
SVector3 normal(double u, double v) const;
void frame(double u, double v, SVector3 &t0, SVector3 &t1, SVector3 &n) const;
void frame(double u, double v,
SPoint3 &p, SVector3 &t0, SVector3 &t1, SVector3 &n) const;
void repositionInnerVertices(const fullMatrix<double> *) const;
};
#endif
......@@ -16,7 +16,7 @@
#include "qualityMeasures.h"
#endif
std::map<int, indicesReversed> MHexahedronN::_order2indicesReversedHex;
std::map<int, IndicesReversed> MHexahedronN::_order2indicesReversedHex;
void MHexahedron::getEdgeRep(bool curved, int num, double *x, double *y, double *z,