Commit 3a255c89 by Amaury Johnen

fix bug reversed hexahedron (issue #293)

parent 5531e2d8
......@@ -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::_order2indicesReversed;
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 _getIndicesReversed(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 = _order2indicesReversed.find(_order);
if (it == _order2indicesReversed.end()) {
indicesReversed indices;
_getIndicesReversed(_order, indices);
_order2indicesReversed[_order] = indices;
it = _order2indicesReversed.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> _order2indicesReversed;
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);
......
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 sign in to comment