Commit 8f498f9f by Christophe Geuzaine

Merge branch 'fixBugsNegativeOrientation' into 'master'

Fix bugs negative orientation

Closes #293

See merge request !24
parents edde29c3 7338fe7b
Pipeline #332 passed with stage
in 9 minutes 0 seconds
......@@ -10,11 +10,14 @@
#include "BasisFactory.h"
#include "polynomialBasis.h"
#include "MQuadrangle.h"
#include "pointsGenerators.h"
#if defined(HAVE_MESH)
#include "qualityMeasures.h"
#endif
std::map<int, indicesReversed> MHexahedronN::_order2indicesReversedHex;
void MHexahedron::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
SVector3 *n)
{
......@@ -407,3 +410,48 @@ int MHexahedronN::getNumFacesRep(bool curved)
return curved ? 6 * (CTX::instance()->mesh.numSubEdges *
CTX::instance()->mesh.numSubEdges * 2) : 12;
}
void _getIndicesReversedHex(int order, indicesReversed &indices)
{
fullMatrix<double> ref = gmshGenerateMonomialsHexahedron(order);
indices.resize(ref.size1());
for (int i = 0; i < ref.size1(); ++i) {
const double u = ref(i, 0);
const double v = ref(i, 1);
const double w = ref(i, 2);
for (int j = 0; j < ref.size1(); ++j) {
if (u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
indices[i] = j;
break;
}
}
}
}
void MHexahedronN::reverse()
{
std::map<int, indicesReversed>::iterator it;
it = _order2indicesReversedHex.find(_order);
if (it == _order2indicesReversedHex.end()) {
indicesReversed indices;
_getIndicesReversedHex(_order, indices);
_order2indicesReversedHex[_order] = indices;
it = _order2indicesReversedHex.find(_order);
}
indicesReversed &indices = it->second;
// copy vertices
std::vector<MVertex*> oldv(8 + _vs.size());
std::copy(_v, _v+8, oldv.begin());
std::copy(_vs.begin(), _vs.end(), oldv.begin()+8);
// reverse
for (int i = 0; i < 8; ++i) {
_v[i] = oldv[indices[i]];
}
for (int i = 0; i < _vs.size(); ++i) {
_vs[i] = oldv[indices[8+i]];
}
}
......@@ -473,8 +473,11 @@ class MHexahedron27 : public MHexahedron {
*
*/
typedef std::vector<int> indicesReversed;
class MHexahedronN : public MHexahedron {
static std::map<int, indicesReversed> _order2indicesReversedHex;
protected:
const char _order;
std::vector<MVertex*> _vs;
......@@ -579,10 +582,7 @@ class MHexahedronN : public MHexahedron {
}
virtual int getNumFacesRep(bool curved);
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
virtual void reverse()
{
Msg::Error("Reverse not implemented yet for MHexahedronN");
}
virtual void reverse();
virtual void getNode(int num, double &u, double &v, double &w) const
{
num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
......
......@@ -7,11 +7,14 @@
#include "Numeric.h"
#include "BasisFactory.h"
#include "Context.h"
#include "pointsGenerators.h"
#if defined(HAVE_MESH)
#include "qualityMeasures.h"
#endif
std::map<int, indicesReversed> MPrismN::_order2indicesReversedPri;
void MPrism::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
SVector3 *n)
{
......@@ -546,3 +549,49 @@ double MPrism::gammaShapeMeasure()
return 0.;
#endif
}
void _getIndicesReversedPri(int order, indicesReversed &indices)
{
fullMatrix<double> ref = gmshGenerateMonomialsPrism(order);
indices.resize(ref.size1());
for (int i = 0; i < ref.size1(); ++i) {
const double u = ref(i, 0);
const double v = ref(i, 1);
const double w = ref(i, 2);
for (int j = 0; j < ref.size1(); ++j) {
if (u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
indices[i] = j;
break;
}
}
}
}
void MPrismN::reverse()
{
std::map<int, indicesReversed>::iterator it;
it = _order2indicesReversedPri.find(_order);
if (it == _order2indicesReversedPri.end()) {
indicesReversed indices;
_getIndicesReversedPri(_order, indices);
_order2indicesReversedPri[_order] = indices;
it = _order2indicesReversedPri.find(_order);
}
indicesReversed &indices = it->second;
// copy vertices
std::vector<MVertex*> oldv(6 + _vs.size());
std::copy(_v, _v+6, oldv.begin());
std::copy(_vs.begin(), _vs.end(), oldv.begin()+6);
// reverse
for (int i = 0; i < 6; ++i) {
_v[i] = oldv[indices[i]];
}
for (int i = 0; i < _vs.size(); ++i) {
_vs[i] = oldv[indices[6+i]];
}
}
......@@ -397,7 +397,12 @@ class MPrism18 : public MPrism {
/*
* MPrismN
*/
typedef std::vector<int> indicesReversed;
class MPrismN : public MPrism {
static std::map<int, indicesReversed> _order2indicesReversedPri;
protected:
std::vector<MVertex *> _vs;
const char _order;
......@@ -504,10 +509,7 @@ class MPrismN : public MPrism {
}
return "";
}
virtual void reverse()
{
Msg::Error("Reverse not implemented yet for MPrismN");
}
virtual void reverse();
virtual void getNode(int num, double &u, double &v, double &w) const
{
const fullMatrix<double> &p = getFunctionSpace()->points;
......
......@@ -8,6 +8,9 @@
#include "Context.h"
#include "BergotBasis.h"
#include "BasisFactory.h"
#include "pointsGenerators.h"
std::map<int, indicesReversed> MPyramidN::_order2indicesReversedPyr;
int MPyramid::getVolumeSign()
{
......@@ -331,3 +334,49 @@ void MPyramidN::getFaceRep(bool curved, int num,
_myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
else MPyramid::getFaceRep(false, num, x, y, z, n);
}
void _getIndicesReversedPyr(int order, indicesReversed &indices)
{
fullMatrix<double> ref = gmshGenerateMonomialsPyramid(order);
indices.resize(ref.size1());
for (int i = 0; i < ref.size1(); ++i) {
const double u = ref(i, 0);
const double v = ref(i, 1);
const double w = ref(i, 2);
for (int j = 0; j < ref.size1(); ++j) {
if (u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
indices[i] = j;
break;
}
}
}
}
void MPyramidN::reverse()
{
std::map<int, indicesReversed>::iterator it;
it = _order2indicesReversedPyr.find(_order);
if (it == _order2indicesReversedPyr.end()) {
indicesReversed indices;
_getIndicesReversedPyr(_order, indices);
_order2indicesReversedPyr[_order] = indices;
it = _order2indicesReversedPyr.find(_order);
}
indicesReversed &indices = it->second;
// copy vertices
std::vector<MVertex*> oldv(5 + _vs.size());
std::copy(_v, _v+5, oldv.begin());
std::copy(_vs.begin(), _vs.end(), oldv.begin()+5);
// reverse
for (int i = 0; i < 5; ++i) {
_v[i] = oldv[indices[i]];
}
for (int i = 0; i < _vs.size(); ++i) {
_vs[i] = oldv[indices[5+i]];
}
}
......@@ -219,7 +219,10 @@ class MPyramid : public MElement {
//------------------------------------------------------------------------------
typedef std::vector<int> indicesReversed;
class MPyramidN : public MPyramid {
static std::map<int, indicesReversed> _order2indicesReversedPyr;
protected:
std::vector<MVertex*> _vs;
......@@ -317,19 +320,7 @@ class MPyramidN : public MPyramid {
Msg::Error("no tag matches a p%d pyramid with %d vertices", _order, 5+_vs.size());
return 0;
}
virtual void reverse()
{
/*
MVertex *tmp;
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
std::vector<MVertex*> inv(_vs.size());
std::vector<int> reverseIndices = _getReverseIndices(_order);
for (unsigned int i = 0; i< _vs.size(); i++)
inv[i] = _vs[reverseIndices[i + 4] - 4];
_vs = inv;
*/
Msg::Error("Reverse not implemented yet for MPyramidN");
}
virtual void reverse();
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z,
SVector3 *n);
virtual int getNumEdgesRep(bool curved);
......
......@@ -350,6 +350,24 @@ double MQuadrangle::getInnerRadius()
#endif // HAVE_LAPACK
}
void MQuadrangleN::reverse()
{
MVertex *tmp;
tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
int npts = _order-1, base = 0;
std::vector<MVertex*>::iterator begin = _vs.begin() + base;
while (npts > 0) {
std::reverse(begin, begin + 4 * npts);
base += 4 * npts;
if (npts > 1) {
tmp = _vs[base+1]; _vs[base+1] = _vs[base+3]; _vs[base+3] = tmp;
}
npts -= 2;
begin = _vs.begin() + base + 4;
}
}
void MQuadrangle::reorient(int rot, bool swap) {
MVertex* tmp[4];
......
......@@ -451,14 +451,7 @@ class MQuadrangleN : public MQuadrangle {
if(_order== 2 && _vs.size() + 4 == 8) return 23;
return MQuadrangle::getTypeForVTK();
}
virtual void reverse()
{
MVertex *tmp;
tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
std::vector<MVertex*> inv;
inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
_vs = inv;
}
virtual void reverse();
// reorient the quadrangle to conform with other face
// orientation computed with MFace based on this face with respect to other
......
......@@ -8,6 +8,7 @@
#include "Numeric.h"
#include "Context.h"
#include "BasisFactory.h"
#include "pointsGenerators.h"
#if defined(HAVE_MESH)
#include "qualityMeasures.h"
......@@ -17,6 +18,8 @@
#define SQU(a) ((a)*(a))
std::map<int, indicesReversed> MTetrahedronN::_order2indicesReversedTet;
void MTetrahedron::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
SVector3 *n)
{
......@@ -315,104 +318,47 @@ void MTetrahedron::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &
Msg::Error("Could not get face information for tetrahedron %d", getNum());
}
static std::vector<std::vector<int> > tetReverseIndices(20);
const std::vector<int> &MTetrahedronN::_getReverseIndices (int order)
void _getIndicesReversedTet(int order, indicesReversed &indices)
{
if(order >= (int)tetReverseIndices.size())
tetReverseIndices.resize(order + 1);
std::vector<int> &r = tetReverseIndices[order];
if (r.size() != 0) return r;
//
// not the funniest code ever ... (guaranteed correct only up to order 5)
//
int nb = (order+1)*(order+2)*(order+3)/6;
r.resize(nb);
int p=0;
for (int layerOrder = order; layerOrder>=0; layerOrder-=4) {
//principal vertices
r[p+0] = p+0;
if (layerOrder ==0) break;
r[p+1] = p+2;
r[p+2] = p+1;
r[p+3] = p+3;
p+=4;
for (int i = 0; i<layerOrder-1; i++) {
//E2 reversed switches with E0
r[p+i] = p+3*(layerOrder-1)-(i+1);
r[p+3*(layerOrder-1)-(i+1)] = p+i;
//E1 is reversed
r[p+(layerOrder-1)+i] = p+2*(layerOrder-1)-(i+1);
//E3 is preserved
r[p+3*(layerOrder-1)+i] = p+3*(layerOrder-1)+i;
//E4 switches with E5
r[p+4*(layerOrder-1)+i] = p+5*(layerOrder-1)+i;
r[p+5*(layerOrder-1)+i] = p+4*(layerOrder-1)+i;
}
p+=6*(layerOrder-1);
//F0(=012) switches its nodes 1 and 2
for (int of = layerOrder-3; of >= 0; of -= 3) {
r[p] = p;
if (of == 0) {
p+=1;
fullMatrix<double> ref = gmshGenerateMonomialsTetrahedron(order);
indices.resize(ref.size1());
for (int i = 0; i < ref.size1(); ++i) {
const double u = ref(i, 0);
const double v = ref(i, 1);
const double w = ref(i, 2);
for (int j = 0; j < ref.size1(); ++j) {
if (u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
indices[i] = j;
break;
}
r[p+1] = p+2;
r[p+2] = p+1;
for (int i = 0; i < of-1; i++) {
//switch edges 0 and 2
r[p+3+i] = p+3+3*(of-1)-(i+1);
r[p+3+3*(of-1)-(i+1)] = p+3+i;
//reverse edge 1
r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1);
}
p += 3*of;
}
//F1 (=013) reversed switches with F2 (=032)
int nf = (layerOrder-2)*(layerOrder-1)/2;
for (int of = layerOrder-3; of >= 0; of -= 3) {
r[p] = p+nf;
r[p+nf] = p;
if (of == 0) {
p += 1;
break;
}
r[p+1] = p+nf+2;
r[p+nf+2] = p+1;
r[p+2] = p+nf+1;
r[p+nf+1] = p+2;
for (int i = 0; i < of-1; i++) {
//switch edges 0 and 2
r[p+3+i] = p+3+3*(of-1)-(i+1)+nf;
r[p+3+3*(of-1)-(i+1)] = p+3+i+nf;
r[p+3+i+nf] = p+3+3*(of-1)-(i+1);
r[p+3+3*(of-1)-(i+1)+nf] = p+3+i;
//reverse edge 1
r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1)+nf;
r[p+3+(of-1)+i+nf] = p+3+2*(of-1)-(i+1);
}
p += 3*of;
}
p+=nf;
//F3(=312) switches its nodes 1 and 2
for (int of = layerOrder-3; of >= 0; of -= 3) {
r[p] = p;
if (of == 0) {
p += 1;
break;
}
r[p+1] = p+2;
r[p+2] = p+1;
for (int i = 0; i < of-1; i++) {
//switch edges 0 and 2
r[p+3+i] = p+3+3*(of-1)-(i+1);
r[p+3+3*(of-1)-(i+1)] = p+3+i;
//reverse edge 1
r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1);
}
p += 3*of;
}
}
return r;
}
void MTetrahedronN::reverse()
{
std::map<int, indicesReversed>::iterator it;
it = _order2indicesReversedTet.find(_order);
if (it == _order2indicesReversedTet.end()) {
indicesReversed indices;
_getIndicesReversedTet(_order, indices);
_order2indicesReversedTet[_order] = indices;
it = _order2indicesReversedTet.find(_order);
}
indicesReversed &indices = it->second;
// copy vertices
std::vector<MVertex*> oldv(4 + _vs.size());
std::copy(_v, _v+4, oldv.begin());
std::copy(_vs.begin(), _vs.end(), oldv.begin()+4);
// reverse
for (int i = 0; i < 4; ++i) {
_v[i] = oldv[indices[i]];
}
for (int i = 0; i < _vs.size(); ++i) {
_vs[i] = oldv[indices[4+i]];
}
}
\ No newline at end of file
......@@ -309,8 +309,11 @@ class MTetrahedron10 : public MTetrahedron {
*
*/
typedef std::vector<int> indicesReversed;
class MTetrahedronN : public MTetrahedron {
static const std::vector<int> &_getReverseIndices(int order);
static std::map<int, indicesReversed> _order2indicesReversedTet;
protected:
std::vector<MVertex *> _vs;
const char _order;
......@@ -412,16 +415,7 @@ class MTetrahedronN : public MTetrahedron {
Msg::Error("no tag matches a p%d tetrahedron with %d vertices", _order, 4+_vs.size());
return 0;
}
virtual void reverse()
{
MVertex *tmp;
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
std::vector<MVertex*> inv(_vs.size());
std::vector<int> reverseIndices = _getReverseIndices(_order);
for (unsigned int i = 0; i< _vs.size(); i++)
inv[i] = _vs[reverseIndices[i + 4] - 4];
_vs = inv;
}
virtual void reverse();
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
virtual int getNumEdgesRep(bool curved);
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
......
......@@ -288,6 +288,25 @@ void MTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*pts = getGQTPts(pOrder);
}
void MTriangleN::reverse()
{
MVertex *tmp;
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
int npts = _order-1, base = 0;
std::vector<MVertex*>::iterator begin = _vs.begin() + base;
while (npts > 0) {
std::reverse(begin, begin + 3 * npts);
base += 3 * npts;
if (npts > 2) {
tmp = _vs[base+1]; _vs[base+1] = _vs[base+2]; _vs[base+2] = tmp;
}
npts -= 3;
begin = _vs.begin() + base + 3;
}
}
void MTriangle::reorient(int rot,bool swap)
{
if (rot == 0 && !swap) return;
......
......@@ -340,14 +340,7 @@ class MTriangleN : public MTriangle {
{
return (_order == 2) ? 22 : MTriangle::getTypeForVTK();
}
virtual void reverse()
{
MVertex *tmp;
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
std::vector<MVertex*> inv;
inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
_vs = inv;
}
virtual void reverse();
virtual void getNode(int num, double &u, double &v, double &w) const
{
num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment