From d63dce3b4ec022742c42382d14caf92539145a66 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Thu, 9 Oct 2014 10:57:04 +0000
Subject: [PATCH] =?UTF-8?q?Add=20FuncSpaceData=20to=20make=20use=20of=20Ba?=
=?UTF-8?q?ses=20simpler=20and=20more=20powerful.=20Add=20a=20validation?=
=?UTF-8?q?=20code=20in=20MetricBasis=20to=20check=20that=20the=20process:?=
=?UTF-8?q?=20"computation=20of=20B=C3=A9zier=20coeff,=20subdivision=20and?=
=?UTF-8?q?=20interpolation"=20is=20correct=20for=20Metric=20(thus=20also?=
=?UTF-8?q?=20for=20Jacobian=20in=203D).?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
Common/GmshDefines.h | 1 +
Geo/MElement.cpp | 141 +++++-
Geo/MElement.h | 18 +-
Geo/MHexahedron.cpp | 1 -
Geo/MHexahedron.h | 6 +-
Geo/MPrism.h | 4 +-
Geo/MPyramid.cpp | 15 +-
Geo/MPyramid.h | 6 +-
Geo/MQuadrangle.h | 2 +-
Geo/MTetrahedron.h | 6 +-
Geo/MTriangle.h | 2 +-
Numeric/BasisFactory.cpp | 65 ++-
Numeric/BasisFactory.h | 62 ++-
Numeric/CMakeLists.txt | 1 +
Numeric/ElementType.cpp | 10 +-
Numeric/FuncSpaceData.cpp | 108 +++++
Numeric/FuncSpaceData.h | 97 +++++
Numeric/JacobianBasis.cpp | 188 +++-----
Numeric/JacobianBasis.h | 49 ++-
Numeric/MetricBasis.cpp | 806 +++++++++++++++++++++++++----------
Numeric/MetricBasis.h | 82 ++--
Numeric/bezierBasis.cpp | 275 +++++++-----
Numeric/bezierBasis.h | 41 +-
Numeric/fullMatrix.h | 10 +-
Numeric/nodalBasis.cpp | 17 +
Numeric/nodalBasis.h | 22 +-
Numeric/pointsGenerators.cpp | 89 +++-
Numeric/pointsGenerators.h | 6 +-
28 files changed, 1504 insertions(+), 626 deletions(-)
create mode 100644 Numeric/FuncSpaceData.cpp
create mode 100644 Numeric/FuncSpaceData.h
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index c85fe0d790..e27349c686 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -66,6 +66,7 @@
#define TYPE_POLYG 9
#define TYPE_POLYH 10
#define TYPE_XFEM 11
+#define TYPE_MINI 12
// Element types in .msh file format (numbers should not be changed)
#define MSH_LIN_2 1
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 82ef462ebe..6a74dd1339 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -49,6 +49,89 @@ MElement::MElement(int num, int part) : _visible(1)
}
}
+MElement* MElement::createElement(int tag, const std::vector<MVertex*> &vertices,
+ int num, int part)
+{
+ const int type = ElementType::ParentTypeFromTag(tag);
+ const int order = ElementType::OrderFromTag(tag);
+ const bool serendipity = ElementType::SerendipityFromTag(tag) > 1;
+
+ if (order == 0) {
+ Msg::Error("p0 elements can not be created (tag %d)", tag);
+ return NULL;
+ }
+
+ switch (type) {
+
+ case TYPE_PNT:
+ return new MPoint(vertices, num, part);
+
+ case TYPE_LIN:
+ if (order == 1)
+ return new MLine(vertices, num, part);
+ else if (order == 2)
+ return new MLine3(vertices, num, part);
+ else
+ return new MLineN(vertices, num, part);
+
+ case TYPE_TRI:
+ if (order == 1)
+ return new MTriangle(vertices, num, part);
+ else if (order == 2)
+ return new MTriangle6(vertices, num, part);
+ else
+ return new MTriangleN(vertices, order, num, part);
+
+ case TYPE_QUA:
+ if (order == 1)
+ return new MQuadrangle(vertices, num, part);
+ else if (order == 2 && serendipity)
+ return new MQuadrangle8(vertices, num, part);
+ else if (order == 2)
+ return new MQuadrangle9(vertices, num, part);
+ else
+ return new MQuadrangleN(vertices, order, num, part);
+
+ case TYPE_TET:
+ if (order == 1)
+ return new MTetrahedron(vertices, num, part);
+ else if (order == 2)
+ return new MTetrahedron10(vertices, num, part);
+ else
+ return new MTetrahedronN(vertices, order, num, part);
+
+ case TYPE_PYR:
+ if (order == 1)
+ return new MPyramid(vertices, num, part);
+ else
+ return new MPyramidN(vertices, order, num, part);
+
+ case TYPE_PRI:
+ if (order == 1)
+ return new MPrism(vertices, num, part);
+ else if (order == 2 && serendipity)
+ return new MPrism15(vertices, num, part);
+ else if (order == 2)
+ return new MPrism18(vertices, num, part);
+ else
+ return new MPrismN(vertices, order, num, part);
+
+ case TYPE_HEX:
+ if (order == 1)
+ return new MHexahedron(vertices, num, part);
+ else if (order == 2 && serendipity)
+ return new MHexahedron20(vertices, num, part);
+ else if (order == 2)
+ return new MHexahedron27(vertices, num, part);
+ else
+ return new MHexahedronN(vertices, order, num, part);
+
+ default:
+ break;
+ }
+ return NULL;
+}
+
void MElement::_getEdgeRep(MVertex *v0, MVertex *v1,
double *x, double *y, double *z, SVector3 *n,
int faceIndex)
@@ -463,8 +546,8 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3]) con
jac[j][1] += ver->y() * gg[j];
jac[j][2] += ver->z() * gg[j];
}
- // printf("GSF (%d,%g %g) = %g %g \n",i,u,v,gg[0],gg[1]);
}
+
return _computeDeterminantAndRegularize(this, jac);
}
@@ -543,6 +626,62 @@ void MElement::getNodesCoord(fullMatrix<double> &nodesXYZ) const
}
}
+double MElement::getEigenvaluesMetric(double u, double v, double w, double values[3]) const
+{
+ double jac[3][3];
+ getJacobian(u, v, w, jac);
+ GradientBasis::mapFromIdealElement(getType(), jac);
+
+ switch (getDim()) {
+ case 1:
+ values[0] = 0;
+ values[1] = -1;
+ values[2] = -1;
+ for (int d = 0; d < 3; ++d)
+ values[0] += jac[d][0] * jac[d][0];
+ return 1;
+
+ case 2:
+ {
+ fullMatrix<double> metric(2, 2);
+ for (int i = 0; i < 2; ++i) {
+ for (int j = 0; j < 2; ++j) {
+ for (int d = 0; d < 3; ++d)
+ metric(i, j) += jac[d][i] * jac[d][j];
+ }
+ }
+ fullVector<double> valReal(values, 2), valImag(2);
+ fullMatrix<double> vecLeft(2, 2), vecRight(2, 2);
+ metric.eig(valReal, valImag, vecLeft, vecRight, true);
+ if (fabs(valImag(0)) + fabs(valImag(1)) > 1e-12)
+ Msg::Warning("not really real");
+ values[2] = -1;
+ return std::sqrt(valReal(0) / valReal(1));
+ }
+
+ case 3:
+ {
+ fullMatrix<double> metric(3, 3);
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ for (int d = 0; d < 3; ++d)
+ metric(i, j) += jac[d][i] * jac[d][j];
+ }
+ }
+
+ fullVector<double> valReal(values, 3), valImag(3);
+ fullMatrix<double> vecLeft(3, 3), vecRight(3, 3);
+ metric.eig(valReal, valImag, vecLeft, vecRight, true);
+
+ return std::sqrt(valReal(0) / valReal(2));
+ }
+
+ default:
+ Msg::Error("wrong dimension for getEigenvaluesMetric function");
+ return -1;
+ }
+}
+
void MElement::pnt(double u, double v, double w, SPoint3 &p) const
{
double x = 0., y = 0., z = 0.;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 26852d821b..5b92db8693 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -18,7 +18,6 @@
#include "polynomialBasis.h"
#include "JacobianBasis.h"
#include "GaussIntegration.h"
-
class GModel;
// A mesh element.
@@ -45,6 +44,10 @@ class MElement
MElement(int num=0, int part=0);
virtual ~MElement(){}
+ // Create an element from tag
+ static MElement* createElement(int tag, const std::vector<MVertex*>&,
+ int num=0, int part=0);
+
// set/get the tolerance for isInside() test
static void setTolerance(const double tol){ _isInsideTolerance = tol; }
static double getTolerance() { return _isInsideTolerance; }
@@ -58,6 +61,15 @@ class MElement
// return the polynomial order the element
virtual int getPolynomialOrder() const { return 1; }
+ // return true if the element can be considered as a serendipity element
+ virtual bool getIsAssimilatedSerendipity() const {
+ return ElementType::SerendipityFromTag(getTypeForMSH()) > 0;
+ }
+ // return true if the element has to be considered as a serendipity element
+ virtual bool getIsOnlySerendipity() const {
+ return ElementType::SerendipityFromTag(getTypeForMSH()) > 1;
+ }
+
// get/set the partition to which the element belongs
virtual int getPartition() const { return _partition; }
virtual void setPartition(int num){ _partition = (short)num; }
@@ -281,6 +293,10 @@ class MElement
virtual const MVertex *getShapeFunctionNode(int i) const { return getVertex(i); }
virtual MVertex *getShapeFunctionNode(int i) { return getVertex(i); }
+ // return the eigenvalues of the metric evaluated at point (u,v,w) in
+ // parametric coordinates
+ virtual double getEigenvaluesMetric(double u, double v, double w, double values[3]) const;
+
// get the point in cartesian coordinates corresponding to the point
// (u,v,w) in parametric coordinates
virtual void pnt(double u, double v, double w, SPoint3 &p) const;
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 39aa47f33d..b8c8c91163 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -7,7 +7,6 @@
#include "Numeric.h"
#include "Context.h"
#include "BasisFactory.h"
-#include "nodalBasis.h"
#include "polynomialBasis.h"
#include "MQuadrangle.h"
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 94fd92b80d..c58701744f 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -486,14 +486,14 @@ class MHexahedronN : public MHexahedron {
virtual int getNumEdgeVertices() const { return 12 * (_order - 1); }
virtual int getNumFaceVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
return 6 * (_order - 1)*(_order - 1);
}
virtual int getNumVolumeVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
return (_order - 1) * (_order - 1) * (_order - 1);
@@ -508,7 +508,7 @@ class MHexahedronN : public MHexahedron {
}
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0) {
+ if (getIsAssimilatedSerendipity()) {
v.resize(4 * _order);
}
else {
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index df37180fcf..b67403cefc 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -413,14 +413,14 @@ class MPrismN : public MPrism {
virtual int getNumEdgeVertices() const { return 9*(_order-1); }
virtual int getNumFaceVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
{int n = _order-1; return (n-1 + 3*n) * n;}
}
virtual int getNumVolumeVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
{int n = _order-1; return n * (n * (n+1) / 2);}
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index 1d2af2da91..968efb6302 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -33,19 +33,12 @@ void MPyramid::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*pts = getGQPyrPts(pOrder);
}
-const JacobianBasis* MPyramid::getJacobianFuncSpace(int o) const
+const JacobianBasis* MPyramid::getJacobianFuncSpace(int order) const
{
- // FIXME add other order and see MPyramid::getFunctionSpace for 'design'
- int order = (o == -1) ? getPolynomialOrder() : o;
-
- switch (order) {
- case 1: return BasisFactory::getJacobianBasis(MSH_PYR_5);
- case 2: return BasisFactory::getJacobianBasis(MSH_PYR_14);
- case 3: return BasisFactory::getJacobianBasis(MSH_PYR_30);
- default: Msg::Error("Order %d pyramid function space not implemented", order); break;
- }
+ if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
- return 0;
+ int tag = ElementType::getTag(TYPE_PYR, order);
+ return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
}
MPyramidN::~MPyramidN() {}
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 7bfdc27caf..b8e7acc708 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -249,7 +249,7 @@ class MPyramidN : public MPyramid {
virtual int getNumEdgeVertices() const { return 8 * (_order - 1); }
virtual int getNumFaceVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
return (_order-1)*(_order-1) + 4 * ((_order - 1) * (_order - 2)) / 2;
@@ -264,7 +264,7 @@ class MPyramidN : public MPyramid {
}
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0) {
+ if (getIsAssimilatedSerendipity()) {
num == 4 ? v.resize(4 * _order)
: v.resize(3 * _order);
}
@@ -290,7 +290,7 @@ class MPyramidN : public MPyramid {
}
virtual int getNumVolumeVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
return (_order-2) * ((_order-2)+1) * (2*(_order-2)+1) / 6;
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 187e6c639f..21b1c29878 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -375,7 +375,7 @@ class MQuadrangleN : public MQuadrangle {
virtual const MVertex *getVertex(int num) const{ return num < 4 ? _v[num] : _vs[num - 4]; }
virtual int getNumFaceVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
return (_order - 1) * (_order - 1);
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 721b88017b..e986c1d037 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -333,7 +333,7 @@ class MTetrahedronN : public MTetrahedron {
virtual int getNumEdgeVertices() const { return 6 * (_order - 1); }
virtual int getNumFaceVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
return 4 * ((_order - 1) * (_order - 2)) / 2;
@@ -348,7 +348,7 @@ class MTetrahedronN : public MTetrahedron {
}
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0) {
+ if (getIsAssimilatedSerendipity()) {
v.resize(3 * _order);
}
else {
@@ -379,7 +379,7 @@ class MTetrahedronN : public MTetrahedron {
}
virtual int getNumVolumeVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
return ((_order - 1) * (_order - 2) * (_order - 3)) / 6;
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 3b8fba8224..83cc1a4d9d 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -278,7 +278,7 @@ class MTriangleN : public MTriangle {
virtual const MVertex *getVertex(int num) const { return num < 3 ? _v[num] : _vs[num - 3]; }
virtual int getNumFaceVertices() const
{
- if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
+ if (getIsAssimilatedSerendipity())
return 0;
else
return (_order - 1) * (_order - 2) / 2;
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index b2906da4c8..d5678f959d 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -3,20 +3,20 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+#include "BasisFactory.h"
#include "GmshDefines.h"
-#include "GmshMessage.h"
-#include "miniBasis.h"
#include "polynomialBasis.h"
#include "pyramidalBasis.h"
-#include "pointsGenerators.h"
-#include "BasisFactory.h"
-#include "MElement.h"
+#include "miniBasis.h"
+#include "MetricBasis.h"
+#include <map>
+#include <cstddef>
std::map<int, nodalBasis*> BasisFactory::fs;
std::map<int, MetricBasis*> BasisFactory::ms;
-BasisFactory::Cont_bezierBasis BasisFactory::bs;
-BasisFactory::Cont_gradBasis BasisFactory::gs;
-BasisFactory::Cont_jacBasis BasisFactory::js;
+std::map<FuncSpaceData, JacobianBasis*> BasisFactory::js;
+std::map<FuncSpaceData, bezierBasis*> BasisFactory::bs;
+std::map<FuncSpaceData, GradientBasis*> BasisFactory::gs;
const nodalBasis* BasisFactory::getNodalBasis(int tag)
{
@@ -67,50 +67,47 @@ const nodalBasis* BasisFactory::getNodalBasis(int tag)
return inserted.first->second;
}
-const JacobianBasis* BasisFactory::getJacobianBasis(int tag, int order)
+const JacobianBasis* BasisFactory::getJacobianBasis(FuncSpaceData fsd)
{
- std::pair<int, int> key(tag, order);
- Cont_jacBasis::iterator it = js.find(key);
- if (it != js.end())
- return it->second;
+ FuncSpaceData data = fsd.getForNonSerendipitySpace();
- JacobianBasis* J = new JacobianBasis(tag, order);
- js.insert(std::make_pair(key, J));
+ std::map<FuncSpaceData, JacobianBasis*>::const_iterator it = js.find(data);
+ if (it != js.end()) return it->second;
+
+ JacobianBasis* J = new JacobianBasis(data);
+ js.insert(std::make_pair(data, J));
return J;
}
const MetricBasis* BasisFactory::getMetricBasis(int tag)
{
std::map<int, MetricBasis*>::const_iterator it = ms.find(tag);
- if (it != ms.end())
- return it->second;
+ if (it != ms.end()) return it->second;
MetricBasis* M = new MetricBasis(tag);
ms.insert(std::make_pair(tag, M));
return M;
}
-const GradientBasis* BasisFactory::getGradientBasis(int tag, int order)
+const GradientBasis* BasisFactory::getGradientBasis(FuncSpaceData data)
{
- std::pair<int, int> key(tag, order);
- Cont_gradBasis::const_iterator it = gs.find(key);
- if (it != gs.end())
- return it->second;
+ std::map<FuncSpaceData, GradientBasis*>::const_iterator it = gs.find(data);
+ if (it != gs.end()) return it->second;
- GradientBasis* G = new GradientBasis(tag, order);
- gs.insert(std::make_pair(key, G));
+ GradientBasis* G = new GradientBasis(data);
+ gs.insert(std::make_pair(data, G));
return G;
}
-const bezierBasis* BasisFactory::getBezierBasis(int parentType, int order)
+const bezierBasis* BasisFactory::getBezierBasis(FuncSpaceData fsd)
{
- std::pair<int, int> key(parentType, order);
- Cont_bezierBasis::iterator it = bs.find(key);
- if (it != bs.end())
- return it->second;
+ FuncSpaceData data = fsd.getForPrimaryElement();
+
+ std::map<FuncSpaceData, bezierBasis*>::const_iterator it = bs.find(data);
+ if (it != bs.end()) return it->second;
- bezierBasis* B = new bezierBasis(parentType, order);
- bs.insert(std::make_pair(key, B));
+ bezierBasis* B = new bezierBasis(data);
+ bs.insert(std::make_pair(data, B));
return B;
}
@@ -130,21 +127,21 @@ void BasisFactory::clearAll()
}
ms.clear();
- Cont_jacBasis::iterator itJ = js.begin();
+ std::map<FuncSpaceData, JacobianBasis*>::iterator itJ = js.begin();
while (itJ != js.end()) {
delete itJ->second;
itJ++;
}
js.clear();
- Cont_gradBasis::iterator itG = gs.begin();
+ std::map<FuncSpaceData, GradientBasis*>::iterator itG = gs.begin();
while (itG != gs.end()) {
delete itG->second;
itG++;
}
gs.clear();
- Cont_bezierBasis::iterator itB = bs.begin();
+ std::map<FuncSpaceData, bezierBasis*>::iterator itB = bs.begin();
while (itB != bs.end()) {
delete itB->second;
itB++;
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index 1125a4b341..9d5ac2acca 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -6,41 +6,65 @@
#ifndef BASISFACTORY_H
#define BASISFACTORY_H
-#include "MElement.h"
-#include "MPyramid.h"
-#include "nodalBasis.h"
-#include "bezierBasis.h"
#include "JacobianBasis.h"
-#include "MetricBasis.h"
+#include "FuncSpaceData.h"
+class nodalBasis;
+class MetricBasis;
+class GradientBasis;
+class bezierBasis;
class BasisFactory
{
- typedef std::map<std::pair<int, int>, JacobianBasis*> Cont_jacBasis;
- typedef std::map<std::pair<int, int>, bezierBasis*> Cont_bezierBasis;
- typedef std::map<std::pair<int, int>, GradientBasis*> Cont_gradBasis;
-
private:
static std::map<int, nodalBasis*> fs;
static std::map<int, MetricBasis*> ms;
- static Cont_jacBasis js;
- static Cont_gradBasis gs;
- static Cont_bezierBasis bs;
- // store bezier bases by parentType and order (no serendipity..)
+ static std::map<FuncSpaceData, JacobianBasis*> js;
+ static std::map<FuncSpaceData, bezierBasis*> bs;
+ static std::map<FuncSpaceData, GradientBasis*> gs;
public:
// Caution: the returned pointer can be NULL
+
+ // Nodal
static const nodalBasis* getNodalBasis(int tag);
- static const JacobianBasis* getJacobianBasis(int tag, int order);
+
+ // Jacobian
+ // Warning: bases returned by BasisFactory::getJacobianBasis(int tag) are the
+ // only safe bases for using Bezier on the jacobian determinant!
+ static const JacobianBasis* getJacobianBasis(FuncSpaceData);
+ static const JacobianBasis* getJacobianBasis(int tag, int order) {
+ const int type = ElementType::ParentTypeFromTag(tag);
+ if (type != TYPE_PYR)
+ return getJacobianBasis(FuncSpaceData(true, tag, order));
+ else
+ return getJacobianBasis(FuncSpaceData(true, tag, false, order+1, order));
+ }
static const JacobianBasis* getJacobianBasis(int tag) {
- return getJacobianBasis(tag, JacobianBasis::jacobianOrder(tag) );
+ const int order = JacobianBasis::jacobianOrder(tag);
+ const int type = ElementType::ParentTypeFromTag(tag);
+ if (type != TYPE_PYR)
+ return getJacobianBasis(FuncSpaceData(true, tag, order));
+ else
+ return getJacobianBasis(FuncSpaceData(true, tag, false, order+2, order));
}
+
+ // Metric
static const MetricBasis* getMetricBasis(int tag);
- static const GradientBasis* getGradientBasis(int tag, int order);
- static const bezierBasis* getBezierBasis(int parentTag, int order);
+ // Gradients
+ static const GradientBasis* getGradientBasis(FuncSpaceData);
+ static const GradientBasis* getGradientBasis(int tag, int order) {
+ return getGradientBasis(FuncSpaceData(true, tag, order));
+ }
+
+ // Bezier
+ static const bezierBasis* getBezierBasis(FuncSpaceData);
+ static const bezierBasis* getBezierBasis(int parentTag, int order) {
+ int primaryTag = ElementType::getTag(parentTag, 1);
+ return getBezierBasis(FuncSpaceData(true, primaryTag, order));
+ }
static const bezierBasis* getBezierBasis(int tag) {
- return getBezierBasis(ElementType::ParentTypeFromTag(tag),
- ElementType::OrderFromTag(tag) );
+ return getBezierBasis(FuncSpaceData(tag));
}
static void clearAll();
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 6977d241ce..a865d9416e 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -7,6 +7,7 @@ set(SRC
Numeric.cpp
fullMatrix.cpp
BasisFactory.cpp
+ FuncSpaceData.cpp
discreteFrechetDistance.cpp
miniBasis.cpp
nodalBasis.cpp
diff --git a/Numeric/ElementType.cpp b/Numeric/ElementType.cpp
index 366a2e1775..7ed6ec3743 100644
--- a/Numeric/ElementType.cpp
+++ b/Numeric/ElementType.cpp
@@ -90,9 +90,11 @@ int ElementType::ParentTypeFromTag(int tag)
case(MSH_PNT_SUB): case(MSH_LIN_SUB):
case(MSH_TRI_SUB): case(MSH_TET_SUB):
return TYPE_XFEM;
+ case(MSH_TRI_MINI):
+ return TYPE_MINI;
default:
- Msg::Error("Unknown element tag %i, assuming tetrahedron.", tag);
- return TYPE_TET;
+ Msg::Error("Unknown element tag %i for parent type, returning -1.", tag);
+ return -1;
}
}
@@ -318,8 +320,8 @@ int ElementType::DimensionFromTag(int tag)
return 3;
default:
- Msg::Error("Unknown element tag %i, assuming dimension 3.", tag);
- return 3;
+ Msg::Error("Unknown element tag %i for dimension, returning -1.", tag);
+ return -1;
}
}
diff --git a/Numeric/FuncSpaceData.cpp b/Numeric/FuncSpaceData.cpp
new file mode 100644
index 0000000000..925443d352
--- /dev/null
+++ b/Numeric/FuncSpaceData.cpp
@@ -0,0 +1,108 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-B-> Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#include "FuncSpaceData.h"
+#include "MElement.h"
+
+FuncSpaceData::FuncSpaceData(MElement *el, int *serendip) :
+ _tag(el->getTypeForMSH()), _spaceOrder(el->getPolynomialOrder()),
+ _nij(0), _nk(_spaceOrder), _pyramidalSpace(el->getType() == TYPE_PYR),
+ _serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
+{}
+
+FuncSpaceData::FuncSpaceData(MElement *el, int order, int *serendip) :
+ _tag(el->getTypeForMSH()), _spaceOrder(order),
+ _nij(0), _nk(_spaceOrder), _pyramidalSpace(el->getType() == TYPE_PYR),
+ _serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
+{}
+
+FuncSpaceData::FuncSpaceData(MElement *el, bool pyr, int nij, int nk, int *serendip) :
+ _tag(el->getTypeForMSH()), _spaceOrder(pyr ? nij+nk : std::max(nij, nk)),
+ _nij(nij), _nk(nk), _pyramidalSpace(pyr),
+ _serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
+{
+ if (el->getType() != TYPE_PYR)
+ Msg::Error("Creation of pyramidal space data for a non-pyramid element !");
+}
+
+FuncSpaceData::FuncSpaceData(int tag, int *serendip) :
+ _tag(tag), _spaceOrder(ElementType::OrderFromTag(tag)),
+ _nij(0), _nk(_spaceOrder),
+ _pyramidalSpace(ElementType::ParentTypeFromTag(tag) == TYPE_PYR),
+ _serendipity(serendip ? *serendip :
+ ElementType::SerendipityFromTag(_tag) > 1)
+{}
+
+FuncSpaceData::FuncSpaceData(bool isTag, int tagOrType, int order,
+ const bool *serendip, bool elemIsSerendip) :
+ _tag(isTag ? tagOrType : ElementType::getTag(tagOrType, order, elemIsSerendip)),
+ _spaceOrder(order), _nij(0), _nk(_spaceOrder),
+ _pyramidalSpace(isTag ?
+ ElementType::ParentTypeFromTag(tagOrType) == TYPE_PYR :
+ tagOrType == TYPE_PYR),
+ _serendipity(serendip ? *serendip :
+ ElementType::SerendipityFromTag(_tag) > 1)
+{}
+
+FuncSpaceData::FuncSpaceData(bool isTag, int tagOrType, bool pyr, int nij,
+ int nk, const bool *serendip, bool elemIsSerendip) :
+ _tag(isTag ? tagOrType :
+ ElementType::getTag(tagOrType, pyr ? nij+nk : std::max(nij, nk), elemIsSerendip)),
+ _spaceOrder(pyr ? nij+nk : std::max(nij, nk)),
+ _nij(nij), _nk(nk), _pyramidalSpace(pyr),
+ _serendipity(serendip ? *serendip :
+ ElementType::SerendipityFromTag(_tag) > 1)
+{
+ if (ElementType::ParentTypeFromTag(_tag) != TYPE_PYR)
+ Msg::Error("Creation of pyramidal space data for a non-pyramid element!");
+}
+
+void FuncSpaceData::getOrderForBezier(int order[3], int exponentZ) const
+{
+ if (_pyramidalSpace && exponentZ < 0) {
+ Msg::Error("getOrderForBezier needs third exponent for pyramidal space!");
+ int a[1];
+ int sum=0;
+ for(int i = 0; i < 1000000; ++i) sum+=a[i];
+ Msg::Info("sum %d", sum);
+ }
+ if (elementType() == TYPE_PYR) {
+ if (_pyramidalSpace) {
+ order[0] = order[1] = _nij + exponentZ;
+ order[2] = _nk;
+ }
+ else {
+ order[0] = order[1] = _nij;
+ order[2] = _nk;
+ }
+ }
+ else
+ order[0] = order[1] = order[2] = _spaceOrder;
+}
+
+FuncSpaceData FuncSpaceData::getForPrimaryElement() const
+{
+ int type = elementType();
+ int primTag = ElementType::getTag(type, 1, elementIsOnlySerendipity());
+
+ if (primTag == _tag) return *this;
+
+ if (type == TYPE_PYR)
+ return FuncSpaceData(true, primTag, _pyramidalSpace, _nij, _nk, &_serendipity);
+ else
+ return FuncSpaceData(true, primTag, _spaceOrder, &_serendipity);
+}
+
+FuncSpaceData FuncSpaceData::getForNonSerendipitySpace() const
+{
+ if (!_serendipity) return *this;
+
+ int type = elementType();
+ bool serendip = false;
+ if (type == TYPE_PYR)
+ return FuncSpaceData(true, _tag, _pyramidalSpace, _nij, _nk, &serendip);
+ else
+ return FuncSpaceData(true, _tag, _spaceOrder, &serendip);
+}
diff --git a/Numeric/FuncSpaceData.h b/Numeric/FuncSpaceData.h
new file mode 100644
index 0000000000..8b6161a76c
--- /dev/null
+++ b/Numeric/FuncSpaceData.h
@@ -0,0 +1,97 @@
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#ifndef FUNCSPACEDATA_H
+#define FUNCSPACEDATA_H
+
+#include "GmshDefines.h"
+#include "GmshMessage.h"
+#include "ElementType.h"
+#include <cstddef>
+class MElement;
+
+class FuncSpaceData
+{
+ // Store data that allow to easily know how to construct gradient, jacobian,
+ // bezier and metric bases.
+
+private:
+ const int _tag, _spaceOrder, _nij, _nk;
+ const bool _pyramidalSpace, _serendipity;
+ // When '_parentType' == TYPE_PYR,
+ // if _pyramidalSpace == true, then the space is {X^i Y^j Z^k | i,j <= k+'_nij', k <= '_nk'},
+ // otherwise, the space is {X^i Y^j Z^k | i,j <= '_nij', k <= '_nk'},
+ // where X = xi/(1-zeta), Y = eta/(1-zeta) and Z = 1-zeta.
+
+public:
+
+ // Constructors using MElement*
+ FuncSpaceData(MElement *el, int *serendip = NULL);
+ FuncSpaceData(MElement *el, int order, int *serendip = NULL);
+ FuncSpaceData(MElement *el, bool pyr, int nij, int nk, int *serendip = NULL);
+
+ // Constructors using element tag
+ FuncSpaceData(int tag, int *serendip = NULL);
+
+ // constructors using element tag or element type
+ FuncSpaceData(bool isTag, int tagOrType, int order,
+ const bool *serendip = NULL, bool elemIsSerendip = false);
+
+ FuncSpaceData(bool isTag, int tagOrType, bool pyr, int nij, int nk,
+ const bool *serendip = NULL, bool elemIsSerendip = false);
+
+ // Print
+ void print() const {
+ Msg::Info("FuncSpaceData: tag%d, order%d, nij%d, nk%d, pyr%d, serendip%d",
+ _tag, _spaceOrder, _nij, _nk, _pyramidalSpace, _serendipity);
+ }
+
+ // Get methods
+ int elementTag() const {return _tag;}
+ int elementType() const {return ElementType::ParentTypeFromTag(_tag);}
+ int elementOrder() const {return ElementType::OrderFromTag(_tag);}
+ int dimension() const {return ElementType::DimensionFromTag(_tag);}
+ int spaceOrder() const {return _spaceOrder;}
+ int nij() const {return _nij;}
+ int nk() const {return _nk;}
+ bool elementIsOnlySerendipity() const {
+ return ElementType::SerendipityFromTag(_tag) > 1;
+ }
+ bool spaceIsSerendipity() const {return _serendipity;}
+ bool isPyramidalSpace() const {return _pyramidalSpace;}
+
+ void getOrderForBezier(int[3], int exponentZ = -1) const;
+
+ // Change space
+ FuncSpaceData getForPrimaryElement() const;
+ FuncSpaceData getForNonSerendipitySpace() const;
+
+ //
+ inline bool operator<(const FuncSpaceData &other) const {
+ if (_tag == other._tag) {
+ if (_spaceOrder == other._spaceOrder) {
+ if (_nij == other._nij) {
+ if (_nk == other._nk) {
+ return _pyramidalSpace == true ? false : other._pyramidalSpace;
+ }
+ else return _nk < other._nk;
+ }
+ else return _nij < other._nij;
+ }
+ else return _spaceOrder < other._spaceOrder;
+ }
+ else return _tag < other._tag;
+ }
+ inline bool operator==(const FuncSpaceData &other) const {
+ return _tag == other._tag && _spaceOrder == other._spaceOrder &&
+ _nij == other._nij && _nk == other._nk &&
+ _pyramidalSpace == other._pyramidalSpace;
+ }
+};
+
+
+
+
+#endif
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index cbfdc0ca27..6a57c63134 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -4,16 +4,11 @@
// bugs and problems to the public mailing list <gmsh@geuz.org>.
#include "JacobianBasis.h"
-
-#include "GmshDefines.h"
-#include "GmshMessage.h"
-
-#include <vector>
-#include "polynomialBasis.h"
-#include "pyramidalBasis.h"
#include "pointsGenerators.h"
+#include "nodalBasis.h"
#include "BasisFactory.h"
#include "Numeric.h"
+#include <cmath>
namespace {
@@ -286,55 +281,24 @@ inline void calcIDI3D(double dxdX, double dxdY, double dxdZ,
}
-GradientBasis::GradientBasis(int tag, int order) :
- _type(ElementType::ParentTypeFromTag(tag))
-{
- const int type = ElementType::ParentTypeFromTag(tag);
+GradientBasis::GradientBasis(FuncSpaceData data)
+ : _data(data) {
fullMatrix<double> samplingPoints;
-
- switch (type) {
- case TYPE_PNT :
- samplingPoints = gmshGeneratePointsLine(0);
- break;
- case TYPE_LIN :
- samplingPoints = gmshGeneratePointsLine(order);
- break;
- case TYPE_TRI :
- samplingPoints = gmshGeneratePointsTriangle(order,false);
- break;
- case TYPE_QUA :
- samplingPoints = gmshGeneratePointsQuadrangle(order,false);
- break;
- case TYPE_TET :
- samplingPoints = gmshGeneratePointsTetrahedron(order,false);
- break;
- case TYPE_PRI :
- samplingPoints = gmshGeneratePointsPrism(order,false);
- break;
- case TYPE_HEX :
- samplingPoints = gmshGeneratePointsHexahedron(order,false);
- break;
- case TYPE_PYR :
- samplingPoints = gmshGeneratePointsPyramidGeneral(false, order+2, order);
- break;
- default :
- Msg::Error("Unknown Jacobian function space for element tag %d", tag);
- return;
- }
+ gmshGeneratePoints(data, samplingPoints);
const int numSampPnts = samplingPoints.size1();
// Store shape function gradients of mapping at Jacobian nodes
fullMatrix<double> allDPsi;
- const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
+ const nodalBasis *mapBasis = BasisFactory::getNodalBasis(_data.elementTag());
mapBasis->df(samplingPoints, allDPsi);
const int numMapNodes = allDPsi.size1();
gradShapeMatX.resize(numSampPnts, numMapNodes);
gradShapeMatY.resize(numSampPnts, numMapNodes);
gradShapeMatZ.resize(numSampPnts, numMapNodes);
- for (int i=0; i<numSampPnts; i++) {
- for (int j=0; j<numMapNodes; j++) {
+ for (int i = 0; i < numSampPnts; i++) {
+ for (int j = 0; j < numMapNodes; j++) {
gradShapeMatX(i, j) = allDPsi(j, 3*i);
gradShapeMatY(i, j) = allDPsi(j, 3*i+1);
gradShapeMatZ(i, j) = allDPsi(j, 3*i+2);
@@ -352,12 +316,13 @@ void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
}
-void GradientBasis::mapFromIdealElement(fullMatrix<double> *dxyzdX,
+void GradientBasis::mapFromIdealElement(int type,
+ fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
- fullMatrix<double> *dxyzdZ) const
+ fullMatrix<double> *dxyzdZ)
{
// 2D scaling
- switch(_type) {
+ switch(type) {
case TYPE_QUA: // Quad, hex, pyramid -> square with side of length 1
case TYPE_HEX:
case TYPE_PYR: {
@@ -374,14 +339,14 @@ void GradientBasis::mapFromIdealElement(fullMatrix<double> *dxyzdX,
}
// 3D scaling
- switch(_type) {
+ switch(type) {
case TYPE_HEX: // Hex, prism -> side of length 1 in z
case TYPE_PRI: {
dxyzdZ->scale(2.);
break;
}
case TYPE_PYR: { // Pyramid -> height sqrt(2)/2
- static const double cPyr = sqrt(2.);
+ static const double cPyr = 1./sqrt(2.);
dxyzdZ->scale(cPyr);
break;
}
@@ -397,55 +362,31 @@ void GradientBasis::mapFromIdealElement(fullMatrix<double> *dxyzdX,
}
}
-JacobianBasis::JacobianBasis(int tag, int jacOrder) :
- _bezier(NULL), _tag(tag), _dim(ElementType::DimensionFromTag(tag)),
- _jacOrder(jacOrder >= 0 ? jacOrder : jacobianOrder(tag))
+void GradientBasis::mapFromIdealElement(int type, double jac[3][3])
{
- const int parentType = ElementType::ParentTypeFromTag(tag);
- const int primJacobianOrder = jacobianOrder(parentType, 1);
+ fullMatrix<double> dxyzdX(jac[0], 1, 3), dxyzdY(jac[1], 1, 3), dxyzdZ(jac[2], 1, 3);
+ mapFromIdealElement(type, &dxyzdX, &dxyzdY, &dxyzdZ);
+}
- fullMatrix<double> lagPoints; // Sampling points
+JacobianBasis::JacobianBasis(FuncSpaceData data)
+ : _bezier(NULL), _data(data), _dim(data.dimension())
+{
+ const int parentType = data.elementType();
+ const int primJacobianOrder = jacobianOrder(parentType, 1);
- switch (parentType) {
- case TYPE_PNT :
- lagPoints = gmshGeneratePointsLine(0);
- break;
- case TYPE_LIN :
- lagPoints = gmshGeneratePointsLine(_jacOrder);
- break;
- case TYPE_TRI :
- lagPoints = gmshGeneratePointsTriangle(_jacOrder,false);
- break;
- case TYPE_QUA :
- lagPoints = gmshGeneratePointsQuadrangle(_jacOrder,false);
- break;
- case TYPE_TET :
- lagPoints = gmshGeneratePointsTetrahedron(_jacOrder,false);
- break;
- case TYPE_PRI :
- lagPoints = gmshGeneratePointsPrism(_jacOrder,false);
- break;
- case TYPE_HEX :
- lagPoints = gmshGeneratePointsHexahedron(_jacOrder,false);
- break;
- case TYPE_PYR :
- lagPoints = gmshGeneratePointsPyramidGeneral(false, _jacOrder+2, _jacOrder);
- break;
- default :
- Msg::Error("Unknown Jacobian function space for element tag %d", tag);
- return;
- }
+ fullMatrix<double> lagPoints; // Sampling points
+ gmshGeneratePoints(data, lagPoints);
numJacNodes = lagPoints.size1();
// Store shape function gradients of mapping at Jacobian nodes
- _gradBasis = BasisFactory::getGradientBasis(tag, _jacOrder);
+ _gradBasis = BasisFactory::getGradientBasis(data);
// Compute matrix for lifting from primary Jacobian basis to Jacobian basis
int primJacType = ElementType::getTag(parentType, primJacobianOrder, false);
const nodalBasis *primJacBasis = BasisFactory::getNodalBasis(primJacType);
numPrimJacNodes = primJacBasis->getNumShapeFunctions();
- matrixPrimJac2Jac.resize(numJacNodes,numPrimJacNodes);
+ matrixPrimJac2Jac.resize(numJacNodes, numPrimJacNodes);
primJacBasis->f(lagPoints, matrixPrimJac2Jac);
// Compute shape function gradients of primary mapping at barycenter,
@@ -471,7 +412,7 @@ JacobianBasis::JacobianBasis(int tag, int jacOrder) :
primGradShapeBarycenterX.resize(numPrimMapNodes);
primGradShapeBarycenterY.resize(numPrimMapNodes);
primGradShapeBarycenterZ.resize(numPrimMapNodes);
- for (int j=0; j<numPrimMapNodes; j++) {
+ for (int j = 0; j < numPrimMapNodes; j++) {
primGradShapeBarycenterX(j) = barDPsi[j][0];
primGradShapeBarycenterY(j) = barDPsi[j][1];
primGradShapeBarycenterZ(j) = barDPsi[j][2];
@@ -480,24 +421,24 @@ JacobianBasis::JacobianBasis(int tag, int jacOrder) :
delete[] barDPsi;
// Compute "fast" Jacobian evaluation matrices (at 1st order nodes + barycenter)
- numJacNodesFast = numPrimMapNodes+1;
- fullMatrix<double> lagPointsFast(numJacNodesFast,3); // Sampling points
- lagPointsFast.copy(primMapBasis->points,0,numPrimMapNodes,
- 0,primMapBasis->points.size2(),0,0); // 1st order nodes
- lagPointsFast(numPrimMapNodes,0) = barycenter[0]; // Last point = barycenter
- lagPointsFast(numPrimMapNodes,1) = barycenter[1];
- lagPointsFast(numPrimMapNodes,2) = barycenter[2];
+ numJacNodesFast = numPrimMapNodes + 1;
+ fullMatrix<double> lagPointsFast(numJacNodesFast, 3); // Sampling points
+ lagPointsFast.copy(primMapBasis->points, 0, numPrimMapNodes,
+ 0, primMapBasis->points.size2(), 0, 0); // 1st order nodes
+ lagPointsFast(numPrimMapNodes, 0) = barycenter[0]; // Last point = barycenter
+ lagPointsFast(numPrimMapNodes, 1) = barycenter[1];
+ lagPointsFast(numPrimMapNodes, 2) = barycenter[2];
fullMatrix<double> allDPsiFast;
- const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
+ const nodalBasis *mapBasis = BasisFactory::getNodalBasis(data.elementTag());
mapBasis->df(lagPointsFast, allDPsiFast);
numMapNodes = mapBasis->getNumShapeFunctions();
gradShapeMatXFast.resize(numJacNodesFast, numMapNodes);
gradShapeMatYFast.resize(numJacNodesFast, numMapNodes);
gradShapeMatZFast.resize(numJacNodesFast, numMapNodes);
- for (int i=0; i<numJacNodesFast; i++) {
- for (int j=0; j<numMapNodes; j++) {
+ for (int i = 0; i < numJacNodesFast; i++) {
+ for (int j = 0; j < numMapNodes; j++) {
gradShapeMatXFast(i, j) = allDPsiFast(j, 3*i);
gradShapeMatYFast(i, j) = allDPsiFast(j, 3*i+1);
gradShapeMatZFast(i, j) = allDPsiFast(j, 3*i+2);
@@ -505,11 +446,11 @@ JacobianBasis::JacobianBasis(int tag, int jacOrder) :
}
}
-const bezierBasis* JacobianBasis::getBezier() const {
- if (_bezier) return _bezier;
- const int parentType = ElementType::ParentTypeFromTag(_tag);
- const_cast<JacobianBasis*>(this)->_bezier =
- BasisFactory::getBezierBasis(parentType, _jacOrder);
+const bezierBasis* JacobianBasis::getBezier() const
+{
+ if (!_bezier) {
+ const_cast<JacobianBasis*>(this)->_bezier = BasisFactory::getBezierBasis(_data);
+ }
return _bezier;
}
@@ -518,10 +459,10 @@ double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullM
{
fullVector<double> dxyzdXbar(3);
- for (int j=0; j<numPrimMapNodes; j++) {
- dxyzdXbar(0) += primGradShapeBarycenterX(j)*nodesXYZ(j,0);
- dxyzdXbar(1) += primGradShapeBarycenterX(j)*nodesXYZ(j,1);
- dxyzdXbar(2) += primGradShapeBarycenterX(j)*nodesXYZ(j,2);
+ for (int j = 0; j < numPrimMapNodes; j++) {
+ dxyzdXbar(0) += primGradShapeBarycenterX(j) * nodesXYZ(j, 0);
+ dxyzdXbar(1) += primGradShapeBarycenterX(j) * nodesXYZ(j, 1);
+ dxyzdXbar(2) += primGradShapeBarycenterX(j) * nodesXYZ(j, 2);
}
if((fabs(dxyzdXbar(0)) >= fabs(dxyzdXbar(1)) && fabs(dxyzdXbar(0)) >= fabs(dxyzdXbar(2))) ||
@@ -667,6 +608,7 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
jacobian.scale(scale);
}
break;
+
}
}
@@ -714,7 +656,7 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
case 0 : {
const int numEl = nodesX.size2();
for (int iEl = 0; iEl < numEl; iEl++)
- for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) = 1.;
+ for (int i = 0; i < nJacNodes; i++) jacobian(i, iEl) = 1.;
break;
}
@@ -860,11 +802,11 @@ inline void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
case 0 : {
for (int i = 0; i < nJacNodes; i++) {
for (int j = 0; j < numMapNodes; j++) {
- JDJ (i,j) = 0.;
- JDJ (i,j+1*numMapNodes) = 0.;
- JDJ (i,j+2*numMapNodes) = 0.;
+ JDJ(i, j) = 0.;
+ JDJ(i, j+1*numMapNodes) = 0.;
+ JDJ(i, j+2*numMapNodes) = 0.;
}
- JDJ(i,3*numMapNodes) = 1.;
+ JDJ(i, 3*numMapNodes) = 1.;
}
break;
}
@@ -1160,9 +1102,9 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
{
// jacobian of the straight elements (only triangles for now)
- SPoint3 v0(nodesXYZ(0,0),nodesXYZ(0,1),nodesXYZ(0,2));
- SPoint3 v1(nodesXYZ(1,0),nodesXYZ(1,1),nodesXYZ(1,2));
- SPoint3 v2(nodesXYZ(2,0),nodesXYZ(2,1),nodesXYZ(2,2));
+ SPoint3 v0(nodesXYZ(0, 0), nodesXYZ(0, 1), nodesXYZ(0, 2));
+ SPoint3 v1(nodesXYZ(1, 0), nodesXYZ(1, 1), nodesXYZ(1, 2));
+ SPoint3 v2(nodesXYZ(2, 0), nodesXYZ(2, 1), nodesXYZ(2, 2));
SPoint3 *IXYZ[3] = {&v0, &v1, &v2};
double jaci[2][2] = {
{IXYZ[1]->x() - IXYZ[0]->x(), IXYZ[2]->x() - IXYZ[0]->x()},
@@ -1174,14 +1116,14 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
for (int l = 0; l < numJacNodes; l++) {
double jac[2][2] = {{0., 0.}, {0., 0.}};
for (int i = 0; i < numMapNodes; i++) {
- const double &dPhidX = _gradBasis->gradShapeMatX(l,i);
- const double &dPhidY = _gradBasis->gradShapeMatY(l,i);
+ const double &dPhidX = _gradBasis->gradShapeMatX(l, i);
+ const double &dPhidY = _gradBasis->gradShapeMatY(l, i);
const double dpsidx = dPhidX * invJaci[0][0] + dPhidY * invJaci[1][0];
const double dpsidy = dPhidX * invJaci[0][1] + dPhidY * invJaci[1][1];
- jac[0][0] += nodesXYZ(i,0) * dpsidx;
- jac[0][1] += nodesXYZ(i,0) * dpsidy;
- jac[1][0] += nodesXYZ(i,1) * dpsidx;
- jac[1][1] += nodesXYZ(i,1) * dpsidy;
+ jac[0][0] += nodesXYZ(i, 0) * dpsidx;
+ jac[0][1] += nodesXYZ(i, 0) * dpsidy;
+ jac[1][0] += nodesXYZ(i, 1) * dpsidx;
+ jac[1][1] += nodesXYZ(i, 1) * dpsidy;
}
const double dxdx = jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1];
const double dydy = jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1];
@@ -1198,8 +1140,8 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
const double aetaxi = ayx * invJaci[0][0] + ayy * invJaci[0][1];
const double axieta = axx * invJaci[1][0] + axy * invJaci[1][1];
for (int i = 0; i < numMapNodes; i++) {
- const double &dPhidX = _gradBasis->gradShapeMatX(l,i);
- const double &dPhidY = _gradBasis->gradShapeMatY(l,i);
+ const double &dPhidX = _gradBasis->gradShapeMatX(l, i);
+ const double &dPhidY = _gradBasis->gradShapeMatY(l, i);
gradLambdaJ(l, i + 0 * numMapNodes) = axixi * dPhidX + axieta * dPhidY;
gradLambdaJ(l, i + 1 * numMapNodes) = aetaxi * dPhidX + aetaeta * dPhidY;
}
@@ -1249,5 +1191,3 @@ int JacobianBasis::jacobianOrder(int parentType, int order)
return 0;
}
}
-
-
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index d45ce88ab8..ace52e9907 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -6,21 +6,19 @@
#ifndef _JACOBIAN_BASIS_H_
#define _JACOBIAN_BASIS_H_
-#include <map>
-#include <vector>
-#include "bezierBasis.h"
#include "fullMatrix.h"
-
+#include "FuncSpaceData.h"
+#include "bezierBasis.h"
class GradientBasis {
- public:
+public:
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
- private:
- const int _type;
+private:
+ const FuncSpaceData _data;
- public:
- GradientBasis(int tag, int order);
+public:
+ GradientBasis(FuncSpaceData);
int getNumSamplingPoints() const {return gradShapeMatX.size1();}
int getNumMapNodes() const {return gradShapeMatX.size2();}
@@ -31,35 +29,42 @@ class GradientBasis {
fullMatrix<double> *dxyzdZ) const;
void mapFromIdealElement(fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
- fullMatrix<double> *dxyzdZ) const;
+ fullMatrix<double> *dxyzdZ) const {
+ GradientBasis::mapFromIdealElement(_data.elementType(), dxyzdX, dxyzdY, dxyzdZ);
+ }
+ static void mapFromIdealElement(int type,
+ fullMatrix<double> *dxyzdX,
+ fullMatrix<double> *dxyzdY,
+ fullMatrix<double> *dxyzdZ);
+ static void mapFromIdealElement(int type, double jac[3][3]);
};
-
class JacobianBasis {
- private:
+private:
const GradientBasis *_gradBasis;
const bezierBasis *_bezier;
- const int _tag, _dim, _jacOrder;
+ const FuncSpaceData _data;
+ const int _dim;
fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
- fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
+ fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
int numJacNodes, numPrimJacNodes;
int numMapNodes, numPrimMapNodes;
int numJacNodesFast;
- public :
- JacobianBasis(int tag, int jacOrder = -1);
+public:
+ JacobianBasis(FuncSpaceData);
// Get methods
- inline int getJacOrder() const { return _jacOrder; }
- inline int getNumJacNodes() const { return numJacNodes; }
- inline int getNumJacNodesFast() const { return numJacNodesFast; }
- inline int getNumMapNodes() const { return numMapNodes; }
- inline int getNumPrimJacNodes() const { return numPrimJacNodes; }
- inline int getNumPrimMapNodes() const { return numPrimMapNodes; }
+ inline int getJacOrder() const {return _data.spaceOrder();}
+ inline int getNumJacNodes() const {return numJacNodes;}
+ inline int getNumJacNodesFast() const {return numJacNodesFast;}
+ inline int getNumMapNodes() const {return numMapNodes;}
+ inline int getNumPrimJacNodes() const {return numPrimJacNodes;}
+ inline int getNumPrimMapNodes() const {return numPrimMapNodes;}
const bezierBasis* getBezier() const;
// Jacobian evaluation methods
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 36cbfbb9db..ce0c248149 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -1,15 +1,16 @@
-// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+#include "MElement.h"
#include "AnalyseCurvedMesh.h"
#include "MetricBasis.h"
#include "BasisFactory.h"
#include "pointsGenerators.h"
#include "BasisFactory.h"
-#include <queue>
#include "OS.h"
+#include <queue>
#include <sstream>
double MetricBasis::_tol = 1e-3;
@@ -48,20 +49,46 @@ namespace {
}
}
-MetricBasis::MetricBasis(int tag)
+MetricBasis::MetricBasis(int tag) : _type(ElementType::ParentTypeFromTag(tag)),
+ _dim(ElementType::DimensionFromTag(tag))
{
- const int type = ElementType::ParentTypeFromTag(tag);
+ const bool serendip = false;
const int metOrder = metricOrder(tag);
- if (type == TYPE_HEX || type == TYPE_PRI) {
- int order = ElementType::OrderFromTag(tag);
- _jacobian = BasisFactory::getJacobianBasis(tag, 3*order);
+ const int jacOrder = 3*metOrder/2;
+
+ // get bezier and gradients for metric space
+ FuncSpaceData *metricSpace = NULL;
+ if (_type != TYPE_PYR)
+ metricSpace = new FuncSpaceData(true, tag, metOrder, &serendip);
+ else
+ metricSpace = new FuncSpaceData(true, tag, false, metOrder+2,
+ metOrder, &serendip);
+
+ _gradients = BasisFactory::getGradientBasis(*metricSpace);
+ _bezier = BasisFactory::getBezierBasis(*metricSpace);
+ delete metricSpace;
+
+ // get jacobian
+ FuncSpaceData *jacSpace = NULL;
+ if (_type == TYPE_TRI || _type == TYPE_QUA) {
+ jacSpace = NULL;
+ }
+ else if (_type == TYPE_TET || _type == TYPE_HEX || _type == TYPE_PRI) {
+ jacSpace = new FuncSpaceData(true, tag, jacOrder, &serendip);
+ }
+ else if (_type == TYPE_PYR) {
+ jacSpace = new FuncSpaceData(true, tag, false, jacOrder+3,
+ jacOrder, &serendip);
}
- else if (type == TYPE_TET || type == TYPE_TRI || type == TYPE_QUA)
- _jacobian = BasisFactory::getJacobianBasis(tag);
else
Msg::Fatal("metric not implemented for element tag %d", tag);
- _gradients = BasisFactory::getGradientBasis(tag, metOrder);
- _bezier = BasisFactory::getBezierBasis(type, metOrder);
+
+ if (jacSpace) {
+ _jacobian = BasisFactory::getJacobianBasis(*jacSpace);
+ delete jacSpace;
+ }
+ else
+ _jacobian = NULL;
_fillInequalities(metOrder);
}
@@ -78,8 +105,18 @@ double MetricBasis::minRCorner(MElement *el)
int order = 1;
if (el->getType() == TYPE_TRI || el->getType() == TYPE_TET) order = 0;
- const GradientBasis *gradients = BasisFactory::getGradientBasis(tag, order);
- const JacobianBasis *jacobian = BasisFactory::getJacobianBasis(tag, order);
+ const GradientBasis *gradients;
+ const JacobianBasis *jacobian;
+ if (el->getType() != TYPE_PYR) {
+ gradients = BasisFactory::getGradientBasis(tag, order);
+ jacobian = BasisFactory::getJacobianBasis(tag, order);
+ }
+ else {
+ const bool serendip = false;
+ FuncSpaceData fsd(true, tag, false, 1, 0, &serendip);
+ gradients = BasisFactory::getGradientBasis(fsd);
+ jacobian = BasisFactory::getJacobianBasis(fsd);
+ }
int nSampPnts = jacobian->getNumJacNodes();
if (el->getType() == TYPE_PYR) nSampPnts = 4;
@@ -144,29 +181,19 @@ double MetricBasis::getMinSampledR(MElement *el, int deg) const
MetricData *md;
_getMetricData(el, md);
- double uvw[3];
- double minmaxQ[2];
- int dim = el->getDim();
- uvw[0] = samplingPoints(0, 0);
- uvw[1] = samplingPoints(0, 1);
- if (dim == 3) uvw[2] = samplingPoints(0, 2);
- else uvw[2] = 0;
- interpolate(el, md, uvw, minmaxQ);
-
- double min, max = min = std::sqrt(minmaxQ[0]/minmaxQ[1]);
- for (int i = 1; i < samplingPoints.size1(); ++i) {
- uvw[0] = samplingPoints(i, 0);
- uvw[1] = samplingPoints(i, 1);
- if (dim == 3) uvw[2] = samplingPoints(i, 2);
- interpolate(el, md, uvw, minmaxQ);
- double tmp = std::sqrt(minmaxQ[0]/minmaxQ[1]);
- min = std::min(min, tmp);
- max = std::max(max, tmp);
- }
+ fullMatrix<double> R;
+ interpolate(el, md, samplingPoints, R);
+
+ if (R.size1() < 1) return -1;
+
+ double min = R(0, 1);
+ for (int i = 1; i < R.size1(); ++i)
+ min = std::min(min, R(i, 1));
+
return min;
}
-double MetricBasis::getBoundMinR(MElement *el)
+double MetricBasis::getBoundMinR(MElement *el) const
{
int nSampPnts = _gradients->getNumSamplingPoints();
int nMapping = _gradients->getNumMapNodes();
@@ -178,12 +205,12 @@ double MetricBasis::getBoundMinR(MElement *el)
// Jacobian coefficients
fullVector<double> jacLag(_jacobian->getNumJacNodes());
fullVector<double> *jac = new fullVector<double>(_jacobian->getNumJacNodes());
- _jacobian->getSignedJacobian(nodes, jacLag);
+ _jacobian->getSignedIdealJacobian(nodes, jacLag);
_jacobian->lag2Bez(jacLag, *jac);
// Metric coefficients
fullMatrix<double> metCoeffLag;
- _fillCoeff<false>(el->getDim(), _gradients, nodes, metCoeffLag);
+ _fillCoeff<true>(el->getDim(), _gradients, nodes, metCoeffLag);
fullMatrix<double> *metCoeff;
metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
_bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
@@ -206,15 +233,21 @@ double MetricBasis::getBoundMinR(MElement *el)
void MetricBasis::_fillInequalities(int metricOrder)
{
- int dimSimplex = _bezier->getDimSimplex();
- int dim = _bezier->getDim();
+ if (_type == TYPE_PYR) {
+ _fillInequalitiesPyr(metricOrder);
+ return;
+ }
+
fullMatrix<int> exp(_bezier->_exponents.size1(), _bezier->_exponents.size2());
for (int i = 0; i < _bezier->_exponents.size1(); ++i) {
for (int j = 0; j < _bezier->_exponents.size2(); ++j) {
exp(i, j) = static_cast<int>(_bezier->_exponents(i, j) + .5);
}
}
+
int ncoeff = _gradients->getNumSamplingPoints();
+ int dimSimplex = _bezier->getDimSimplex();
+ int dim = _bezier->getDim();
int countP3 = 0, countJ2 = 0, countA = 0;
for (int i = 0; i < ncoeff; i++) {
@@ -296,47 +329,160 @@ void MetricBasis::_fillInequalities(int metricOrder)
}
}
- if (_jacobian) {
- exp.resize(_jacobian->getBezier()->_exponents.size1(),
- _jacobian->getBezier()->_exponents.size2());
- for (int i = 0; i < _jacobian->getBezier()->_exponents.size1(); ++i) {
- for (int j = 0; j < _jacobian->getBezier()->_exponents.size2(); ++j) {
- exp(i, j) = static_cast<int>(_jacobian->getBezier()->_exponents(i, j) + .5);
+ if (_dim == 2) {
+ _lightenInequalities(countJ2, countP3, countA);
+ return;
+ }
+
+ exp.resize(_jacobian->getBezier()->_exponents.size1(),
+ _jacobian->getBezier()->_exponents.size2());
+ for (int i = 0; i < _jacobian->getBezier()->_exponents.size1(); ++i) {
+ for (int j = 0; j < _jacobian->getBezier()->_exponents.size2(); ++j) {
+ exp(i, j) = static_cast<int>(_jacobian->getBezier()->_exponents(i, j) + .5);
+ }
+ }
+ int njac = _jacobian->getNumJacNodes();
+ for (int i = 0; i < njac; i++) {
+ for (int j = i; j < njac; j++) {
+ int order = metricOrder/2*3;
+ double num = 1, den = 1;
+ {
+ int compl1 = order;
+ int compl2 = order;
+ int compltot = 2*order;
+ for (int k = 0; k < dimSimplex; k++) {
+ num *= nChoosek(compl1, exp(i, k))
+ * nChoosek(compl2, exp(j, k));
+ den *= nChoosek(compltot, exp(i, k) + exp(j, k));
+ compl1 -= exp(i, k);
+ compl2 -= exp(j, k);
+ compltot -= exp(i, k) + exp(j, k);
+ }
}
+ for (int k = dimSimplex; k < dim; k++) {
+ num *= nChoosek(order, exp(i, k))
+ * nChoosek(order, exp(j, k));
+ den *= nChoosek(2*order, exp(i, k) + exp(j, k));
+ }
+
+ if (i != j) num *= 2;
+
+ ++countJ2;
+ int hash = 0;
+ for (int k = 0; k < dim; k++) {
+ hash += (exp(i, k)+exp(j, k)) * pow_int(2*order+1, k);
+ }
+ _ineqJ2[hash].push_back(IneqData(num/den, i, j));
}
- int njac = _jacobian->getNumJacNodes();
- for (int i = 0; i < njac; i++) {
- for (int j = i; j < njac; j++) {
- int order = metricOrder/2*3;
- double num = 1, den = 1;
- {
- int compl1 = order;
- int compl2 = order;
- int compltot = 2*order;
- for (int k = 0; k < dimSimplex; k++) {
- num *= nChoosek(compl1, exp(i, k))
- * nChoosek(compl2, exp(j, k));
- den *= nChoosek(compltot, exp(i, k) + exp(j, k));
- compl1 -= exp(i, k);
- compl2 -= exp(j, k);
- compltot -= exp(i, k) + exp(j, k);
- }
+ }
+
+ _lightenInequalities(countJ2, countP3, countA);
+}
+
+void MetricBasis::_fillInequalitiesPyr(int metricOrder)
+{
+ fullMatrix<int> exp(_bezier->_exponents.size1(), _bezier->_exponents.size2());
+ for (int i = 0; i < _bezier->_exponents.size1(); ++i) {
+ for (int j = 0; j < _bezier->_exponents.size2(); ++j) {
+ exp(i, j) = static_cast<int>(_bezier->_exponents(i, j) + .5);
+ }
+ }
+
+ int ncoeff = _gradients->getNumSamplingPoints();
+
+ int countP3 = 0, countJ2 = 0, countA = 0;
+ for (int i = 0; i < ncoeff; i++) {
+ for (int j = i; j < ncoeff; j++) {
+ double num, den;
+ num = nChoosek(metricOrder+2, exp(i, 0))
+ * nChoosek(metricOrder+2, exp(i, 1))
+ * nChoosek(metricOrder , exp(i, 2))
+ * nChoosek(metricOrder+2, exp(j, 0))
+ * nChoosek(metricOrder+2, exp(j, 1))
+ * nChoosek(metricOrder , exp(j, 2));
+ den = nChoosek(2*metricOrder+4, exp(i, 0) + exp(j, 0))
+ * nChoosek(2*metricOrder+4, exp(i, 1) + exp(j, 1))
+ * nChoosek(2*metricOrder , exp(i, 2) + exp(j, 2));
+
+ if (i != j) num *= 2;
+
+ ++countA;
+ int hash = 0;
+ for (int l = 0; l < 3; l++) {
+ hash += (exp(i, l)+exp(j, l)) * pow_int(2*metricOrder+1, l);
+ }
+ _ineqA[hash].push_back(IneqData(num/den, i, j));
+
+
+ for (int k = j; k < ncoeff; ++k) {
+ double num, den;
+ num = nChoosek(metricOrder+2, exp(i, 0))
+ * nChoosek(metricOrder+2, exp(i, 1))
+ * nChoosek(metricOrder , exp(i, 2))
+ * nChoosek(metricOrder+2, exp(j, 0))
+ * nChoosek(metricOrder+2, exp(j, 1))
+ * nChoosek(metricOrder , exp(j, 2))
+ * nChoosek(metricOrder+2, exp(k, 0))
+ * nChoosek(metricOrder+2, exp(k, 1))
+ * nChoosek(metricOrder , exp(k, 2));
+ den = nChoosek(3*metricOrder+6, exp(i, 0) + exp(j, 0) + exp(k, 0))
+ * nChoosek(3*metricOrder+6, exp(i, 1) + exp(j, 1) + exp(k, 1))
+ * nChoosek(3*metricOrder , exp(i, 2) + exp(j, 2) + exp(k, 2));
+
+ if (i == j) {
+ if (j != k) num *= 3;
}
- for (int k = dimSimplex; k < dim; k++) {
- num *= nChoosek(order, exp(i, k))
- * nChoosek(order, exp(j, k));
- den *= nChoosek(2*order, exp(i, k) + exp(j, k));
+ else {
+ if (j == k || i == k) {
+ num *= 3;
+ }
+ else num *= 6;
}
- if (i != j) num *= 2;
-
- ++countJ2;
+ ++countP3;
int hash = 0;
- for (int k = 0; k < dim; k++) {
- hash += (exp(i, k)+exp(j, k)) * pow_int(2*order+1, k);
+ for (int l = 0; l < 3; l++) {
+ hash += (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*metricOrder+1, l);
}
- _ineqJ2[hash].push_back(IneqData(num/den, i, j));
+ if (j == k && j != i)
+ _ineqP3[hash].push_back(IneqData(num/den, k, j, i));
+ else
+ _ineqP3[hash].push_back(IneqData(num/den, i, j, k));
+ }
+ }
+ }
+
+ exp.resize(_jacobian->getBezier()->_exponents.size1(),
+ _jacobian->getBezier()->_exponents.size2());
+ for (int i = 0; i < _jacobian->getBezier()->_exponents.size1(); ++i) {
+ for (int j = 0; j < _jacobian->getBezier()->_exponents.size2(); ++j) {
+ exp(i, j) = static_cast<int>(_jacobian->getBezier()->_exponents(i, j) + .5);
+ }
+ }
+ int njac = _jacobian->getNumJacNodes();
+ for (int i = 0; i < njac; i++) {
+ for (int j = i; j < njac; j++) {
+ int order = metricOrder/2*3;
+
+ double num, den;
+ num = nChoosek(order+3, exp(i, 0))
+ * nChoosek(order+3, exp(i, 1))
+ * nChoosek(order , exp(i, 2))
+ * nChoosek(order+3, exp(j, 0))
+ * nChoosek(order+3, exp(j, 1))
+ * nChoosek(order , exp(j, 2));
+ den = nChoosek(2*order+6, exp(i, 0) + exp(j, 0))
+ * nChoosek(2*order+6, exp(i, 1) + exp(j, 1))
+ * nChoosek(2*order , exp(i, 2) + exp(j, 2));
+
+ if (i != j) num *= 2;
+
+ ++countJ2;
+ int hash = 0;
+ for (int k = 0; k < 3; k++) {
+ hash += (exp(i, k)+exp(j, k)) * pow_int(2*order+1, k);
}
+ _ineqJ2[hash].push_back(IneqData(num/den, i, j));
}
}
@@ -378,150 +524,297 @@ void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
counta -= cnt[2];
}
-void MetricBasis::interpolate(const MElement *el, const MetricData *md, const double *uvw, double *minmaxQ) const
+namespace {
+ static double symRand(double f = 1)
+ {
+ return f * (rand()%2001 - 1000) / 1000.;
+ }
+}
+
+bool MetricBasis::validateBezierForMetricAndJacobian()
{
- if (minmaxQ == NULL) {
- Msg::Error("Cannot write solution of interpolation");
- return;
+ Msg::Info("Testing Bezier interpolation and subdivision "
+ "for jacobien and metric on all element types...");
+ int numError = 0;
+
+ // Parameters
+ const int numType = 6;
+ const int acceptedTypes[numType] = {TYPE_TRI, TYPE_QUA, TYPE_TET,
+ TYPE_PYR, TYPE_PRI, TYPE_HEX};
+ const int maxOrder = 3; // at least 3 (so that serendip tet are tested)
+ const int numElem = 5; // at least 2 (first is reference element)
+ const int numSubdiv = 2; // at least 1
+ const int numSampPnt = 10; // at least 1
+ const double toleranceTensor = 1e-11; // at most 1e-5 (metric takes values in [0,1])
+ double tolerance; // computed in function of tag
+
+ //
+ static const double epsilon = std::numeric_limits<double>::epsilon();
+ double sumRatio = 0;
+ int numRatio = 0;
+
+ for (int tag = 1; tag <= MSH_NUM_TYPE; ++tag) {
+ if (tag > 66 && tag < 71) continue; //not conventional elements
+ if (tag > 75 && tag < 79) continue; //no element tag 76, 77, 78...
+
+ // Check if accepted type
+ const int type = ElementType::ParentTypeFromTag(tag);
+ bool knownType = false;
+ for (int i = 0; i < numType; ++i) {
+ knownType = (knownType || type == acceptedTypes[i]);
+ }
+ if (!knownType) continue;
+
+ const int order = ElementType::OrderFromTag(tag);
+ const int dim = ElementType::DimensionFromTag(tag);
+ const bool serendip = ElementType::SerendipityFromTag(tag) > 1;
+
+ // Skip p0 elements and elements for which order > 'maxOrder'
+ // and skip for now serendipty pyramids (not implemented)
+ if (order < 1 || order > maxOrder) continue;
+ if (type == TYPE_PYR && serendip) continue;
+
+ Msg::Info("... testing element tag %d", tag);
+
+ // Compute tolerance
+ tolerance = epsilon * pow_int(10, order*dim);
+ if (type == TYPE_PYR) tolerance = std::max(tolerance, epsilon*1e9);
+
+ // Get reference nodes
+ const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
+ fullMatrix<double> nodes;
+ mapBasis->getReferenceNodes(nodes);
+
+ // Create 'numElem' elements more and more randomized
+ for (int iel = 0; iel < numElem; ++iel) {
+ const double range = static_cast<double>(iel) / (numElem-1) / order;
+ std::vector<MVertex*> vertices(nodes.size1());
+ for (int i = 0; i < nodes.size1(); ++i) {
+ vertices[i] = new MVertex(nodes(i, 0) + symRand(range),
+ dim > 1 ? nodes(i, 1) + symRand(range) : 0,
+ dim > 2 ? nodes(i, 2) + symRand(range) : 0);
+ }
+ MElement *el = MElement::createElement(tag, vertices);
+ if (!el) {
+ Msg::Error("MElement was unable to create element for tag %d", tag);
+ ++numError;
+ }
+
+ const MetricBasis *metricBasis = BasisFactory::getMetricBasis(tag);
+
+ // compare the two metric
+
+ fullMatrix<double> metric_Bez(numSampPnt, 2);
+ fullVector<int> isub(numSubdiv);
+ fullMatrix<double> uvw(numSampPnt, 3);
+ metricBasis->interpolateAfterNSubdivisions(el, numSubdiv, numSampPnt,
+ isub, uvw, metric_Bez);
+
+ double sumTol = 0;
+ int numBadMatch = 0;
+ int numBadMatchTensor = 0;
+ double maxBadMatch = 0;
+ double maxBadMatchTensor = 0;
+ for (int isamp = 0; isamp < numSampPnt; ++isamp) {
+ double dum[3];
+ double &u = uvw(isamp, 0);
+ double &v = uvw(isamp, 1);
+ double &w = uvw(isamp, 2);
+ double metric_Lag = el->getEigenvaluesMetric(u, v, w, dum);
+
+ double diff = std::abs(metric_Lag - metric_Bez(isamp, 0));
+ double diffTensor = std::abs(metric_Lag - metric_Bez(isamp, 1));
+ sumTol += diff;
+ double ratio = (metric_Bez(isamp, 0)-metric_Lag) / tolerance;
+ sumRatio += ratio;
+
+ if (diffTensor > toleranceTensor) {
+ ++numBadMatchTensor;
+ maxBadMatchTensor = std::max(maxBadMatchTensor, diffTensor);
+ }
+ else if (diff > tolerance) {
+ ++numBadMatch;
+ maxBadMatch = std::max(maxBadMatch, diff);
+ }
+ }
+ if (numBadMatchTensor > .2*numSampPnt) {
+ ++numError;
+ Msg::Error("Too much errors "
+ "even when computing by metric tensor (max %g)", maxBadMatchTensor);
+ }
+ else if (numBadMatch > .5*numSampPnt) {
+ ++numError;
+ Msg::Error("Too much errors (max %g)", maxBadMatch);
+ }
+ }
}
- int order = _bezier->getOrder();
+ if (numError) Msg::Error("Validation of Bezier terminated with %d errors, "
+ "you have work...", numError);
+ else Msg::Info("Validation of Bezier terminated without errors", numError);
+ return numError;
+}
- int dimSimplex = 0;
- fullMatrix<double> exponents;
- double bezuvw[3];
- switch (el->getType()) {
- case TYPE_PYR:
- bezuvw[0] = .5 * (1 + uvw[0]);
- bezuvw[1] = .5 * (1 + uvw[1]);
- bezuvw[2] = uvw[2];
- //_interpolateBezierPyramid(uvw, minmaxQ);
- return;
+void MetricBasis::interpolate(const MElement *el,
+ const MetricData *md,
+ const fullMatrix<double> &nodes,
+ fullMatrix<double> &R) const
+{
+ fullMatrix<double> &metcoeffs = *md->_metcoeffs, *metric = new fullMatrix<double>;
+ fullVector<double> &jaccoeffs = *md->_jaccoeffs, *jac = new fullVector<double>;
- case TYPE_HEX:
- bezuvw[0] = .5 * (1 + uvw[0]);
- bezuvw[1] = .5 * (1 + uvw[1]);
- bezuvw[2] = .5 * (1 + uvw[2]);
- dimSimplex = 0;
- exponents = gmshGenerateMonomialsHexahedron(order);
- break;
+ _bezier->interpolate(metcoeffs, nodes, *metric);
- case TYPE_TET:
- bezuvw[0] = uvw[0];
- bezuvw[1] = uvw[1];
- bezuvw[2] = uvw[2];
- dimSimplex = 3;
- exponents = gmshGenerateMonomialsTetrahedron(order);
+ R.resize(nodes.size1(), 2);
+
+ switch (metcoeffs.size2()) {
+ case 1:
+ for (int i = 0; i < R.size1(); ++i)
+ R(i, 0) = R(i, 1) = 1;
break;
- case TYPE_PRI:
- bezuvw[0] = uvw[0];
- bezuvw[1] = uvw[1];
- bezuvw[2] = .5 * (1 + uvw[2]);
- dimSimplex = 2;
- exponents = gmshGenerateMonomialsPrism(order);
+ case 3:
+ for (int i = 0; i < R.size1(); ++i) {
+ // Compute from q, p
+ double p = pow((*metric)(i, 1), 2);
+ p += pow((*metric)(i, 2), 2);
+ p = std::sqrt(p);
+ R(i, 0) = std::sqrt(_R2Dsafe((*metric)(i, 0), p));
+ // Comppute from tensor
+ fullMatrix<double> metricTensor(2, 2);
+ metricTensor(0, 0) = (*metric)(i, 0) + (*metric)(i, 1);
+ metricTensor(1, 1) = (*metric)(i, 0) - (*metric)(i, 1);
+ metricTensor(0, 1) = metricTensor(1, 0) = (*metric)(i, 2);
+ fullVector<double> valReal(2), valImag(2);
+ fullMatrix<double> vecLeft(2, 2), vecRight(2, 2);
+ metricTensor.eig(valReal, valImag, vecLeft, vecRight, true);
+ R(i, 1) = std::sqrt(valReal(0) / valReal(1));
+ }
break;
- case TYPE_TRI:
- bezuvw[0] = uvw[0];
- bezuvw[1] = uvw[1];
- bezuvw[2] = uvw[2];
- dimSimplex = 2;
- exponents = gmshGenerateMonomialsTriangle(order);
+ case 7:
+ _jacobian->getBezier()->interpolate(jaccoeffs, nodes, *jac);
+ for (int i = 0; i < R.size1(); ++i) {
+ // Compute from q, p, J
+ double p = 0;
+ for (int k = 1; k < 7; ++k) p += pow((*metric)(i, k), 2);
+ p = std::sqrt(p);
+ R(i, 0) = std::sqrt(_R3Dsafe((*metric)(i, 0), p, (*jac)(i)));
+ // Comppute from tensor
+ fullMatrix<double> metricTensor(3, 3);
+ for (int k = 0; k < 3; ++k) {
+ static double fact1 = std::sqrt(6);
+ static double fact2 = std::sqrt(3);
+ const int ki = k%2;
+ const int kj = std::min(k+1, 2);
+ metricTensor(k, k) = (*metric)(i, k+1)*fact1 + (*metric)(i, 0);
+ metricTensor(ki, kj) = metricTensor(kj, ki) = (*metric)(i, k+4)*fact2;
+ }
+ fullVector<double> valReal(3), valImag(3);
+ fullMatrix<double> vecLeft(3, 3), vecRight(3, 3);
+ metricTensor.eig(valReal, valImag, vecLeft, vecRight, true);
+ R(i, 1) = std::sqrt(valReal(0) / valReal(2));
+ }
break;
+
+ default:
+ Msg::Error("Wrong number of functions for metric: %d",
+ metcoeffs.size2());
}
- int numCoeff = exponents.size1();
- int dim = exponents.size2();
-
- fullMatrix<double> metcoeffs = *md->_metcoeffs;
- fullVector<double> jaccoeffs = *md->_jaccoeffs;
-
- double *terms = new double[metcoeffs.size2()];
- for (int t = 0; t < metcoeffs.size2(); ++t) {
- terms[t] = 0;
- for (int i = 0; i < numCoeff; i++) {
- double dd = 1;
- double pointCompl = 1.;
- int exponentCompl = order;
- for (int k = 0; k < dimSimplex; k++) {
- dd *= nChoosek(exponentCompl, (int) exponents(i, k))
- * pow(bezuvw[k], exponents(i, k));
- pointCompl -= bezuvw[k];
- exponentCompl -= (int) exponents(i, k);
- }
- dd *= pow(pointCompl, exponentCompl);
+ delete jac;
+ //delete metric;
+}
- for (int k = dimSimplex; k < dim; k++)
- dd *= nChoosek(order, (int) exponents(i, k))
- * pow(bezuvw[k], exponents(i, k))
- * pow(1. - bezuvw[k], order - exponents(i, k));
- terms[t] += metcoeffs(i, t) * dd;
- }
+void MetricBasis::interpolateAfterNSubdivisions(
+ const MElement *el, int numSubdiv, int numPnt,
+ fullVector<int> &isub,
+ fullMatrix<double> &uvw,
+ fullMatrix<double> &metric) const
+{
+ // Interpolate metric after 'numSub' random subdivision at
+ // 'numPnt' random points.
+ // Return: isub, the subdomain tag of each subdivision,
+ // uvw, the reference points at which metric has been interpolated.
+ // metric, the interpolation.
+
+ MetricData *md;
+ _getMetricData(el, md);
+
+ // Keep trace of subdomain to be able to compute uvw.
+ // (Ensure to have the tag for the complete element):
+ const nodalBasis *mapBasis = BasisFactory::getNodalBasis(el->getTypeForMSH());
+ fullMatrix<double> nodes;
+ mapBasis->getReferenceNodesForBezier(nodes);
+
+ const int nSub = _bezier->getNumDivision();
+ const int numCoeff = md->_metcoeffs->size2();
+ const int numMetPnts = md->_metcoeffs->size1();
+ const int numJacPnts = md->_jaccoeffs ? md->_jaccoeffs->size() : 0;
+ const int numNodPnts = nodes.size1();
+
+ const bezierBasis *bezierMapping;
+ if (el->getType() != TYPE_PYR)
+ bezierMapping = BasisFactory::getBezierBasis(el->getTypeForMSH());
+ else {
+ FuncSpaceData data(true, el->getTypeForMSH(), false,
+ el->getPolynomialOrder(), el->getPolynomialOrder());
+ bezierMapping = BasisFactory::getBezierBasis(data);
}
- switch (metcoeffs.size2()) {
- case 1:
- minmaxQ[0] = terms[0];
- minmaxQ[1] = terms[0];
- break;
+ for (int k = 0; k < numSubdiv; ++k) {
+ fullMatrix<double> subcoeffs, subnodes;
+ fullVector<double> subjac;
+ _bezier->subdivideBezCoeff(*md->_metcoeffs, subcoeffs);
+ bezierMapping->subdivideBezCoeff(nodes, subnodes);
- case 3:
- {
- double tmp = pow(terms[1], 2);
- tmp += pow(terms[2], 2);
- tmp = std::sqrt(tmp);
- minmaxQ[0] = terms[0] - tmp;
- minmaxQ[1] = terms[0] + tmp;
+ if (_dim == 3)
+ _jacobian->getBezier()->subdivideBezCoeff(*md->_jaccoeffs, subjac);
+ delete md;
+
+ isub(k) = std::rand() % nSub;
+
+ fullMatrix<double> *coeff = new fullMatrix<double>(numMetPnts, numCoeff);
+ coeff->copy(subcoeffs, isub(k) * numMetPnts, numMetPnts, 0, numCoeff, 0, 0);
+ nodes.copy(subnodes, isub(k) * numNodPnts, numNodPnts, 0, _dim, 0, 0);
+ fullVector<double> *jac = NULL;
+ if (_dim == 3) {
+ jac = new fullVector<double>(numJacPnts);
+ jac->copy(subjac, isub(k) * numJacPnts, numJacPnts, 0);
+ }
+
+ md = new MetricData(coeff, jac);
}
- break;
- case 7:
+ // compute a random convex combination of reference nodes
+ fullMatrix<double> subuvw(uvw.size1(), _dim);
{
- double tmp = pow(terms[1], 2);
- tmp += pow(terms[2], 2);
- tmp += pow(terms[3], 2);
- tmp += pow(terms[4], 2);
- tmp += pow(terms[5], 2);
- tmp += pow(terms[6], 2);
- tmp = std::sqrt(tmp);
- double factor = std::sqrt(6)/3;
- if (tmp < 1e-3*terms[0]) {
- minmaxQ[0] = terms[0] - factor * tmp;
- minmaxQ[1] = terms[0] + factor * tmp;
- }
- else {
- double phi;
- //{
- fullMatrix<double> nodes(1, 3);
- nodes(0, 0) = uvw[0];
- nodes(0, 1) = uvw[1];
- nodes(0, 2) = uvw[2];
-
- fullMatrix<double> result;
- _jacobian->interpolate(jaccoeffs, nodes, result, true);
- phi = result(0, 0)*result(0, 0);
- //}
- phi -= terms[0]*terms[0]*terms[0];
- phi += .5*terms[0]*tmp*tmp;
- phi /= tmp*tmp*tmp;
- phi *= 3*std::sqrt(6);
- if (phi > 1) phi = 1;
- if (phi < -1) phi = -1;
- phi = std::acos(phi)/3;
- minmaxQ[0] = terms[0] + factor * tmp * std::cos(phi + 2*M_PI/3);
- minmaxQ[1] = terms[0] + factor * tmp * std::cos(phi);
- ((MetricBasis*)this)->file << terms[0] << " " << tmp/std::sqrt(6) << " " << result(0, 0) << std::endl;
+ int tagPrimary = ElementType::getTag(el->getType(), 1);
+ const nodalBasis *primMapBasis = BasisFactory::getNodalBasis(tagPrimary);
+ fullMatrix<double> refNodes;
+ primMapBasis->getReferenceNodes(refNodes);
+ double *c = new double[refNodes.size1()];
+
+ for (int k = 0; k < uvw.size1(); ++k) {
+ double sum = 0;
+ int exp = 1 + std::rand() % 5;
+ for (int i = 0; i < refNodes.size1(); ++i) {
+ c[i] = pow_int((std::rand() % 1000) / 1000., exp);
+ sum += c[i];
+ }
+ for (int i = 0; i < refNodes.size1(); ++i) {
+ c[i] /= sum;
+ subuvw(k, 0) += c[i] * refNodes(i, 0);
+ if (_dim > 1) subuvw(k, 1) += c[i] * refNodes(i, 1);
+ if (_dim > 2) subuvw(k, 2) += c[i] * refNodes(i, 2);
+ }
}
- }
- break;
- default:
- Msg::Error("Wrong number of functions for metric: %d",
- metcoeffs.size2());
+ delete[] c;
}
- delete[] terms;
+ interpolate(el, md, subuvw, metric);
+ bezierMapping->interpolate(nodes, subuvw, uvw, false);
}
int MetricBasis::metricOrder(int tag)
@@ -535,12 +828,13 @@ int MetricBasis::metricOrder(int tag)
case TYPE_LIN : return order;
case TYPE_TRI :
- case TYPE_TET : return 2*order-2;
+ case TYPE_TET :
+ case TYPE_PYR : return 2*order-2;
case TYPE_QUA :
case TYPE_PRI :
- case TYPE_HEX :
- case TYPE_PYR : return 2*order;
+ case TYPE_HEX : return 2*order;
+
default :
Msg::Error("Unknown element type %d, return order 0", parentType);
return 0;
@@ -688,7 +982,7 @@ void MetricBasis::_computeRmax(
for (int k = 1; k < 7; ++k) {
p += pow_int(coeff(i, k), 2);
}
- p = std::sqrt(p/6);
+ p = std::sqrt(p);
const double a = q/p;
if (a > 1e4) {
RmaxLag = std::max(RmaxLag, std::sqrt((a - std::sqrt(3)) / (a + std::sqrt(3))));
@@ -791,7 +1085,7 @@ void MetricBasis::_computeTermBeta(double &a, double &K,
dRda = sin * sqrt + .5 * term1 * (1-a*a);
}
-void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
+void MetricBasis::_getMetricData(const MElement *el, MetricData *&md) const
{
int nSampPnts = _gradients->getNumSamplingPoints();
int nMapping = _gradients->getNumMapNodes();
@@ -801,14 +1095,17 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
el->getNodesCoord(nodes);
// Jacobian coefficients
- fullVector<double> jacLag(_jacobian->getNumJacNodes());
- fullVector<double> *jac = new fullVector<double>(_jacobian->getNumJacNodes());
- _jacobian->getSignedJacobian(nodes, jacLag);
- _jacobian->lag2Bez(jacLag, *jac);
+ fullVector<double> *jac = NULL;
+ if (_dim == 3) {
+ fullVector<double> jacLag(_jacobian->getNumJacNodes());
+ jac = new fullVector<double>(_jacobian->getNumJacNodes());
+ _jacobian->getSignedIdealJacobian(nodes, jacLag);
+ _jacobian->lag2Bez(jacLag, *jac);
+ }
// Metric coefficients
fullMatrix<double> metCoeffLag;
- _fillCoeff<false>(el->getDim(), _gradients, nodes, metCoeffLag);
+ _fillCoeff<true>(el->getDim(), _gradients, nodes, metCoeffLag);
fullMatrix<double> *metCoeff;
metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
_bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
@@ -818,7 +1115,7 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
template<bool ideal>
void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
- fullMatrix<double> &nodes, fullMatrix<double> &coeff)
+ const fullMatrix<double> &nodes, fullMatrix<double> &coeff)
{
const int nSampPnts = gradients->getNumSamplingPoints();
@@ -837,8 +1134,17 @@ void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
coeff.resize(nSampPnts, 3);
for (int i = 0; i < nSampPnts; i++) {
- const double &dxdX = dxydX(i,0), &dydX = dxydX(i,1), &dzdX = dxydX(i,2);
- const double &dxdY = dxydY(i,0), &dydY = dxydY(i,1), &dzdY = dxydY(i,2);
+ const double &dxdX = dxydX(i,0), &dydX = dxydX(i,1);
+ const double &dxdY = dxydY(i,0), &dydY = dxydY(i,1);
+ double dzdX, dzdY;
+ if (nodes.size2() > 2) {
+ dzdX = dxydX(i,2);
+ dzdY = dxydY(i,2);
+ }
+ else {
+ dzdX = 0;
+ dzdY = 0;
+ }
const double dvxdX = dxdX*dxdX + dydX*dydX + dzdX*dzdX;
const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY;
coeff(i, 0) = (dvxdX + dvxdY) / 2;
@@ -863,13 +1169,14 @@ void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY;
const double dvxdZ = dxdZ*dxdZ + dydZ*dydZ + dzdZ*dzdZ;
coeff(i, 0) = (dvxdX + dvxdY + dvxdZ) / 3;
- coeff(i, 1) = dvxdX - coeff(i, 0);
- coeff(i, 2) = dvxdY - coeff(i, 0);
- coeff(i, 3) = dvxdZ - coeff(i, 0);
- const double fact = std::sqrt(2);
- coeff(i, 4) = fact * (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
- coeff(i, 5) = fact * (dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY);
- coeff(i, 6) = fact * (dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ);
+ static double fact1 = 1./std::sqrt(6);
+ static double fact2 = 1./std::sqrt(3);
+ coeff(i, 1) = fact1 * (dvxdX - coeff(i, 0));
+ coeff(i, 2) = fact1 * (dvxdY - coeff(i, 0));
+ coeff(i, 3) = fact1 * (dvxdZ - coeff(i, 0));
+ coeff(i, 4) = fact2 * (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
+ coeff(i, 5) = fact2 * (dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY);
+ coeff(i, 6) = fact2 * (dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ);
}
}
break;
@@ -891,7 +1198,7 @@ double MetricBasis::_computeMinlagR(const fullVector<double> &jac,
for (int k = 1; k < 7; ++k) {
p += pow_int(coeff(i, k), 2);
}
- p = std::sqrt(p/6);
+ p = std::sqrt(p);
const double a = q/p;
if (a > 1e4) { // TODO: from _tol ?
Rmin = std::min(Rmin, std::sqrt((a - std::sqrt(3)) / (a + std::sqrt(3))));
@@ -907,10 +1214,9 @@ double MetricBasis::_computeMinlagR(const fullVector<double> &jac,
for (int i = 0; i < num; ++i) {
if (jac(i) <= 0.) return 0;
- const double q = coeff(i, 0);
+ const double &q = coeff(i, 0);
const double p = pow_int(coeff(i, 1), 2) + pow_int(coeff(i, 2), 2);
- const double a = q/std::sqrt(p);
- const double tmpR = _R2Dsafe(a);
+ const double tmpR = _R2Dsafe(q, std::sqrt(p));
Rmin = std::min(Rmin, std::sqrt(tmpR));
}
return Rmin;
@@ -945,10 +1251,6 @@ void MetricBasis::_minMaxA(
max = std::max(val, max);
++it;
}
- if (coeff.size2() == 7) {
- min *= 6;
- max *= 6;
- }
min = min > 1 ? std::sqrt(min) : 1;
max = std::sqrt(max);
@@ -963,7 +1265,7 @@ void MetricBasis::_minK(const fullMatrix<double> &coeff,
for (int l = 1; l < 7; ++l) {
r(i) += coeff(i, l) * coeff(i, l);
}
- r(i) = std::sqrt(r(i)/6);
+ r(i) = std::sqrt(r(i));
}
min = 1e10;
@@ -1015,7 +1317,7 @@ void MetricBasis::_maxAstKpos(const fullMatrix<double> &coeff,
for (int l = 1; l < 7; ++l) {
P(i) += coeff(i, l) * coeff(i, l);
}
- P(i) = std::sqrt(P(i)/6);
+ P(i) = std::sqrt(P(i));
}
double min = 1e10;
@@ -1057,13 +1359,12 @@ void MetricBasis::_maxAstKneg(const fullMatrix<double> &coeff,
for (int l = 1; l < 7; ++l) {
P(i) += coeff(i, l) * coeff(i, l);
}
- P(i) = std::sqrt(P(i)/6);
+ P(i) = std::sqrt(P(i));
for (int j = 0; j < coeff.size1(); ++j) {
Q(i, j) = 0;
for (int l = 1; l < 7; ++l) {
Q(i, j) += coeff(i, l) * coeff(j, l);
}
- Q(i, j) /= 6;
}
}
@@ -1108,7 +1409,7 @@ void MetricBasis::_maxKstAfast(const fullMatrix<double> &coeff,
for (int l = 1; l < 7; ++l) {
r(i) += coeff(i, l) * coeff(i, l);
}
- r(i) = std::sqrt(r(i)/6);
+ r(i) = std::sqrt(r(i));
}
double min = 1e10;
@@ -1150,13 +1451,12 @@ void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff,
for (int l = 1; l < 7; ++l) {
P(i) += coeff(i, l) * coeff(i, l);
}
- P(i) = std::sqrt(P(i)/6);
+ P(i) = std::sqrt(P(i));
for (int j = 0; j < coeff.size1(); ++j) {
Q(i, j) = 0;
for (int l = 1; l < 7; ++l) {
Q(i, j) += coeff(i, l) * coeff(j, l);
}
- Q(i, j) /= 6;
}
}
@@ -1191,3 +1491,59 @@ void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff,
maxK = 1/beta*(mina*mina*mina-min);
}
+
+double MetricBasis::_R3Dsafe(double q, double p, double J)
+{
+ if (q > 1e5*p) {
+ const double m = p*std::sqrt(3)/q;
+ return (1-m) / (1+m);
+ }
+ const double a = q/p;
+ const double K = J*J/p/p/p;
+ //Msg::Info("%g %g", a-3, K-16);
+ return _R3Dsafe(a, K);
+}
+
+double MetricBasis::_R3Dsafe(double a, double K)
+{
+ const double x = .5 * (K + (3 - a*a)*a);
+ if (x > 1+1e-7 || x < -1-1e-7) {
+ Msg::Warning("x = %g (a,K) = (%g,%g)", x, a, K);
+ //Msg::Fatal("a");
+ }
+
+ double ans;
+ if (x >= 1) ans = (a - 1) / (a + 2);
+ else if (x <= -1) ans = (a - 2) / (a + 1);
+ else {
+ const double phi = std::acos(x) / 3;
+ //Msg::Warning("phi %g", phi - M_PI/3);
+ ans = (a + 2*std::cos(phi + 2*M_PI/3)) / (a + 2*std::cos(phi));
+ }
+
+ if (ans < 0 || ans > 1) {
+ //Msg::Warning("R = %g", ans);
+ if (ans < 0) return 0;
+ else return 1;
+ }
+ return ans;
+}
+
+double MetricBasis::_R2Dsafe(double q, double p)
+{
+ if (q < 0 || p < 0 || p > q) {
+ Msg::Error("wrong argument for 2d metric (%g, %g)", q, p);
+ int a[1];
+ int sum=0;
+ for(int i = 0; i < 1000000; ++i) sum+=a[i];
+ Msg::Info("sum %d", sum);
+ }
+ return (q-p) / (q+p);
+}
+
+double MetricBasis::_R2Dsafe(double a)
+{
+ if (a < 1 || !std::isfinite(a))
+ Msg::Error("wrong argument for 2d metric (%g)", a);
+ return (a - 1) / (a + 1);
+}
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index fc15f45a5c..2416f8d12f 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -6,20 +6,22 @@
#ifndef _METRIC_BASIS_H_
#define _METRIC_BASIS_H_
-#include "MElement.h"
#include "JacobianBasis.h"
#include "fullMatrix.h"
#include <fstream>
#include <cmath>
+class MElement;
class MetricBasis {
friend class MetricCoefficient;
friend class GMSH_AnalyseCurvedMeshPlugin;
+
private:
const JacobianBasis *_jacobian;
const GradientBasis *_gradients;
const bezierBasis *_bezier;
static double _tol;
+ const int _type, _dim;
std::fstream file;
@@ -32,6 +34,8 @@ private:
IneqData(double val, int i, int j, int k = -1) : i(i), j(j), k(k), val(val) {}
};
+ std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3, _ineqA;
+
class MetricData {
public:
fullMatrix<double> *_metcoeffs;
@@ -40,7 +44,8 @@ private:
int _depth, _num;
public:
- MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d, int num) :
+ MetricData(fullMatrix<double> *m, fullVector<double> *j,
+ double r = -1, int d = -1, int num = -1) :
_metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d), _num(num) {}
~MetricData() {
delete _metcoeffs;
@@ -48,25 +53,50 @@ private:
}
};
- std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3, _ineqA;
-
public:
MetricBasis(int elementTag);
static void setTol(double tol) {_tol = tol;}
+
+ const JacobianBasis* getJacobianForMetric() const {return _jacobian;}
+ const bezierBasis* getBezier() const {return _bezier;}
+ int getNumMetricNodes() const {return _gradients->getNumSamplingPoints();}
+
static double boundMinR(MElement *el);
- static double minRCorner(MElement *el);
static double minSampledR(MElement *el, int order);
+ double getBoundMinR(MElement*) const;
+ double getMinSampledR(MElement*, int order) const;
- double getBoundMinR(MElement*);
- double getMinSampledR(MElement*, int) const;
+ static double minRCorner(MElement *el);
- void interpolate(const MElement*, const MetricData*, const double *uvw, double *minmaxQ) const;
+ template<bool ideal>
+ void getMetricCoeff(const fullMatrix<double> &nodes,
+ fullMatrix<double> &coeff) const {
+ _fillCoeff<ideal>(_dim, _gradients, nodes, coeff);
+ }
+ void lag2Bez(const fullMatrix<double> &metCoeffLag,
+ fullMatrix<double> &metCoeffBez) const {
+ _bezier->matrixLag2Bez.mult(metCoeffLag, metCoeffBez);
+ }
static int metricOrder(int tag);
+public:
+ // Validation for computation of Bezier coefficients & subdivision
+ // of Jacobian determinant and Metric stuffs
+ static bool validateBezierForMetricAndJacobian();
+ void interpolate(const MElement*, const MetricData*, const double *uvw, double *minmaxQ) const;
+ void interpolate(const MElement*, const MetricData*,
+ const fullMatrix<double> &nodes, fullMatrix<double> &R) const;
+ void interpolateAfterNSubdivisions(const MElement *el,
+ int numSubdiv, int numPnt,
+ fullVector<int> &isub,
+ fullMatrix<double> &uvw,
+ fullMatrix<double> &metric) const;
+
private:
void _fillInequalities(int order);
+ void _fillInequalitiesPyr(int order);
void _lightenInequalities(int&, int&, int&); //TODO change
void _computeRmin(const fullMatrix<double>&, const fullVector<double>&,
@@ -75,12 +105,12 @@ private:
double &RmaxLag) const;
void _computeTermBeta(double &a, double &K, double &dRda,
double &term1, double &phip) const;
- void _getMetricData(MElement*, MetricData*&) const;
+ void _getMetricData(const MElement*, MetricData*&) const;
double _subdivideForRmin(MetricData*, double RminLag, double tol) const;
template<bool ideal>
static void _fillCoeff(int dim, const GradientBasis*,
- fullMatrix<double> &nodes, fullMatrix<double> &coeff);
+ const fullMatrix<double> &nodes, fullMatrix<double> &coeff);
static double _computeMinlagR(const fullVector<double> &jac,
const fullMatrix<double> &coeff, int num);
@@ -95,32 +125,10 @@ private:
void _maxKstAsharp(const fullMatrix<double>&, const fullVector<double>&,
double mina, double beta, double &maxK) const;
- static double _R3Dsafe(double a, double K) {
- const double x = .5 * (K - a*a*a + 3*a);
- if (x > 1+1e-13 || x < -1-1e-13) {
- Msg::Warning("x = %g (|1+%g|)", x, std::abs(x)-1);
- Msg::Fatal("a");
- }
-
- double ans;
- if (x >= 1) ans = (a - 1) / (a + 2);
- else if (x <= -1) ans = (a - 2) / (a + 1);
- else {
- const double phi = std::acos(x) / 3;
- ans = (a + 2*std::cos(phi + 2*M_PI/3)) / (a + 2*std::cos(phi));
- }
-
- if (ans < 0 || ans > 1) {
- Msg::Warning("R = %g", ans);
- if (ans < 0) return 0;
- else return 1;
- }
- return ans;
- }
- static double _R2Dsafe(double a) {
- if (a < 1) Msg::Warning("R2d = %g", (a - 1) / (a + 1));
- return (a - 1) / (a + 1);
- }
+ static double _R3Dsafe(double a, double K);
+ static double _R3Dsafe(double q, double p, double J);
+ static double _R2Dsafe(double a);
+ static double _R2Dsafe(double q, double p);
private:
class gterIneq {
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index fc8a48f060..5a372bfdb7 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -12,6 +12,7 @@
#include "pointsGenerators.h"
#include "BasisFactory.h"
#include "Numeric.h"
+#include <sstream>
namespace {
@@ -224,14 +225,14 @@ std::vector< fullMatrix<double> > generateSubPointsHex(int order)
return subPoints;
}
-std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
+std::vector< fullMatrix<double> > generateSubPointsPyr(int nij, int nk)
{
- if (order == 0) {
+ if (nk == 0) {
std::vector< fullMatrix<double> > subPoints(4);
fullMatrix<double> prox;
- subPoints[0] = gmshGenerateMonomialsPyramidGeneral(false, 2, 0);
- subPoints[0].scale(.5/2);
+ subPoints[0] = gmshGenerateMonomialsPyramidGeneral(false, nij, nk);
+ subPoints[0].scale(.5/nij);
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -251,11 +252,11 @@ std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
std::vector< fullMatrix<double> > subPoints(8);
fullMatrix<double> ref, prox;
- subPoints[0] = gmshGenerateMonomialsPyramidGeneral(false, order+2, order);
+ subPoints[0] = gmshGenerateMonomialsPyramidGeneral(false, nij, nk);
prox.setAsProxy(subPoints[0], 2, 1);
prox.scale(-1);
- prox.add(order);
- subPoints[0].scale(.5/(order+2));
+ prox.add(nk);
+ subPoints[0].scale(.5/std::max(nij,nk));
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -285,13 +286,10 @@ std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
prox.setAsProxy(subPoints[7], 2, 1);
prox.add(.5);
- const int nPts = subPoints[0].size1();
for (int i = 0; i < 8; ++i) {
- for (int j = 0; j < nPts; ++j) {
- const double factor = (1. - subPoints[i](j, 2));
- subPoints[i](j, 0) = subPoints[i](j, 0) * factor;
- subPoints[i](j, 1) = subPoints[i](j, 1) * factor;
- }
+ prox.setAsProxy(subPoints[i], 2, 1);
+ prox.scale(-1);
+ prox.add(1);
}
return subPoints;
@@ -303,6 +301,10 @@ int nChoosek(int n, int k)
{
if (n < k || k < 0) {
Msg::Error("Wrong argument for combination. (%d, %d)", n, k);
+ int a[10];
+ int e = 0;
+ for (int i = 0; i < 1000000; ++i) e += a[i];
+ Msg::Info("e%d", e);
return 1;
}
@@ -361,33 +363,35 @@ fullMatrix<double> generateBez2LagMatrix
fullMatrix<double> generateBez2LagMatrixPyramid
(const fullMatrix<double> &exponent, const fullMatrix<double> &point,
- int order)
+ bool pyr, int nij, int nk)
{
if(exponent.size1() != point.size1() || exponent.size2() != point.size2() ||
exponent.size2() != 3){
Msg::Fatal("Wrong sizes for pyramid's bez2lag matrix generation %d %d -- %d %d",
- exponent.size1(),point.size1(),
- exponent.size2(),point.size2());
+ exponent.size1(), point.size1(),
+ exponent.size2(), point.size2());
return fullMatrix<double>(1, 1);
}
const int ndofs = exponent.size1();
+ int n01 = nij;
fullMatrix<double> bez2Lag(ndofs, ndofs);
for (int i = 0; i < ndofs; i++) {
for (int j = 0; j < ndofs; j++) {
- double denom = 1. - point(i, 2);
+ if (pyr) n01 = exponent(j, 2) + nij;
bez2Lag(i, j) =
- nChoosek(order + 2, exponent(j, 0))
- * nChoosek(order + 2, exponent(j, 1))
- * nChoosek(order, exponent(j, 2))
- * pow_int(point(i, 0) / denom, exponent(j, 0))
- * pow_int(point(i, 1) / denom, exponent(j, 1))
- * pow_int(1. - point(i, 2) , exponent(j, 2))
- * pow_int(1. - point(i, 0) / denom, order + 2 - exponent(j, 0))
- * pow_int(1. - point(i, 1) / denom, order + 2 - exponent(j, 1))
- * pow_int(point(i, 2) , order - exponent(j, 2));
+ nChoosek(n01, exponent(j, 0))
+ * nChoosek(n01, exponent(j, 1))
+ * nChoosek(nk , exponent(j, 2))
+ * pow_int(point(i, 0), exponent(j, 0))
+ * pow_int(point(i, 1), exponent(j, 1))
+ * pow_int(point(i, 2), exponent(j, 2))
+ * pow_int(1. - point(i, 0), n01 - exponent(j, 0))
+ * pow_int(1. - point(i, 1), n01 - exponent(j, 1))
+ * pow_int(1. - point(i, 2), nk - exponent(j, 2));
}
+ //Msg::Info("bb %g %g %g", point(i, 0) / denom, point(i, 1) / denom, 1. - point(i, 2));
}
return bez2Lag;
}
@@ -420,7 +424,7 @@ fullMatrix<double> generateSubDivisor
fullMatrix<double> generateSubDivisorPyramid
(const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
- const fullMatrix<double> &lag2Bez, int order)
+ const fullMatrix<double> &lag2Bez, bool pyr, int nij, int nk)
{
if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
@@ -437,7 +441,8 @@ fullMatrix<double> generateSubDivisorPyramid
for (unsigned int i = 0; i < subPoints.size(); i++) {
fullMatrix<double> intermediate1 =
- generateBez2LagMatrixPyramid(exponents, subPoints[i], order);
+ generateBez2LagMatrixPyramid(exponents, subPoints[i],
+ pyr, nij, nk);
lag2Bez.mult(intermediate1, intermediate2);
subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
}
@@ -445,56 +450,79 @@ fullMatrix<double> generateSubDivisorPyramid
}
}
-void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
- const fullMatrix<double> &uvw,
- fullMatrix<double> &result,
- bool bezCoord) const
+void bezierBasis::generateBezierPoints(fullMatrix<double> &points) const
{
- result.resize(uvw.size1(), coeffs.size2());
- int dimSimplex;
- fullMatrix<double> bezuvw = uvw;
- switch (_type) {
- case TYPE_HEX:
- if (!bezCoord) {
- bezuvw.add(1);
- bezuvw.scale(.5);
- }
- dimSimplex = 0;
- break;
+ gmshGenerateMonomials(_data, points);
+ points.scale(1./_data.spaceOrder());
+
+ if (_data.elementType() == TYPE_PYR && _data.nk() < _data.spaceOrder()) {
+ fullMatrix<double> prox;
+ prox.setAsProxy(points, 2, 1);
+ prox.add(1-static_cast<double>(_data.nk())/_data.spaceOrder());
+ }
+}
+void bezierBasis::_FEpoints2BezPoints(fullMatrix<double> &points) const
+{
+ switch (_data.elementType()) {
+ case TYPE_TRI:
case TYPE_TET:
- dimSimplex = 3;
+ break;
+
+ case TYPE_QUA:
+ case TYPE_HEX:
+ points.add(1);
+ points.scale(.5);
break;
case TYPE_PRI:
- if (!bezCoord) {
+ {
fullMatrix<double> tmp;
- tmp.setAsProxy(bezuvw, 3, 1);
+ tmp.setAsProxy(points, 2, 1);
tmp.add(1);
tmp.scale(.5);
}
- dimSimplex = 2;
break;
- default:
case TYPE_PYR:
- Msg::Error("Bezier interpolation not implemented for type of element %d", _type);
- /*bezuvw[0] = .5 * (1 + uvw[0]);
- bezuvw[1] = .5 * (1 + uvw[1]);
- bezuvw[2] = uvw[2];
- _interpolateBezierPyramid(uvw, minmaxQ);*/
+ for (int i = 0; i < points.size1(); ++i) {
+ points(i, 2) = 1. - points(i, 2);
+ points(i, 0) = .5 * (1 + points(i, 0) / points(i, 2));
+ points(i, 1) = .5 * (1 + points(i, 1) / points(i, 2));
+ }
+ break;
+
+ default:
+ Msg::Error("_FEpoints2BezPoints not implemented for "
+ "type of element %d", _data.elementType());
return;
}
+}
- int numCoeff = _exponents.size1();
+void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
+ const fullMatrix<double> &uvw,
+ fullMatrix<double> &result,
+ bool bezCoord) const
+{
+ if (result.size1() < uvw.size1() || result.size2() < coeffs.size2())
+ result.resize(uvw.size1(), coeffs.size2());
+
+ fullMatrix<double> bezuvw = uvw;
+ if (!bezCoord) _FEpoints2BezPoints(bezuvw);
+
+ const int numCoeff = _exponents.size1();
+ const int dim = _exponents.size2();
+ int order[3];
for (int m = 0; m < uvw.size1(); ++m) {
for (int n = 0; n < coeffs.size2(); ++n) result(m, n) = 0;
for (int i = 0; i < numCoeff; i++) {
+ _data.getOrderForBezier(order, _exponents(i, dim-1));
+ //if (m==0) Msg::Info("i=%d exp %g %g %g order %d %d %d", i, _exponents(i, 0), _exponents(i, 1), _exponents(i, 2), order[0], order[1], order[2]);
double dd = 1;
double pointCompl = 1.;
- int exponentCompl = _order;
- for (int k = 0; k < dimSimplex; k++) {
+ int exponentCompl = order[0];
+ for (int k = 0; k < _dimSimplex; k++) {
dd *= nChoosek(exponentCompl, (int) _exponents(i, k))
* pow(bezuvw(m, k), _exponents(i, k));
pointCompl -= bezuvw(m, k);
@@ -502,10 +530,10 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
}
dd *= pow_int(pointCompl, exponentCompl);
- for (int k = dimSimplex; k < _exponents.size2(); k++) {
- dd *= nChoosek(_order, (int) _exponents(i, k))
+ for (int k = _dimSimplex; k < dim; k++) {
+ dd *= nChoosek(order[k], (int) _exponents(i, k))
* pow_int(bezuvw(m, k), _exponents(i, k))
- * pow_int(1. - bezuvw(m, k), _order - _exponents(i, k));
+ * pow_int(1. - bezuvw(m, k), order[k] - _exponents(i, k));
}
for (int n = 0; n < coeffs.size2(); ++n)
result(m, n) += coeffs(i, n) * dd;
@@ -513,36 +541,35 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
}
}
-void bezierBasis::_construct(int parentType, int p)
+void bezierBasis::subdivideBezCoeff(const fullMatrix<double> &coeff,
+ fullMatrix<double> &subCoeff) const
{
- _order = p;
- _type = parentType;
- std::vector< fullMatrix<double> > subPoints;
-
- if (parentType == TYPE_PYR) {
- _numLagCoeff = p == 0 ? 4 : 8;
- _exponents = gmshGenerateMonomialsPyramidGeneral(false, _order+2, _order);
-
- subPoints = generateSubPointsPyr(_order);
- _numDivisions = static_cast<int>(subPoints.size());
+ if (subCoeff.size1() < subDivisor.size1()
+ || subCoeff.size2() < coeff.size2() ) {
+ subCoeff.resize(subDivisor.size1(), coeff.size2());
+ }
+ subDivisor.mult(coeff, subCoeff);
+}
- fullMatrix<double> bezierPoints;
- bezierPoints.resize(_exponents.size1(), _exponents.size2());
- const double ord = _order + 2;
- for (int i = 0; i < bezierPoints.size1(); ++i) {
- bezierPoints(i, 2) = (_order - _exponents(i, 2)) / ord;
- const double scale = 1. - bezierPoints(i, 2);
- bezierPoints(i, 0) = _exponents(i, 0) / ord * scale;
- bezierPoints(i, 1) = _exponents(i, 1) / ord * scale;
- }
+void bezierBasis::subdivideBezCoeff(const fullVector<double> &coeff,
+ fullVector<double> &subCoeff) const
+{
+ if (subCoeff.size() < subDivisor.size1()) {
+ subCoeff.resize(subDivisor.size1());
+ }
+ subDivisor.mult(coeff, subCoeff);
+}
- matrixBez2Lag = generateBez2LagMatrixPyramid(_exponents, bezierPoints, _order);
- matrixBez2Lag.invert(matrixLag2Bez);
- subDivisor = generateSubDivisorPyramid(_exponents, subPoints, matrixLag2Bez, _order);
- return;
+void bezierBasis::_construct()
+{
+ if (_data.elementType() == TYPE_PYR) {
+ Msg::Fatal("This bezierBasis constructor is not for pyramids !");
}
- switch (parentType) {
+ std::vector< fullMatrix<double> > subPoints;
+ int order = _data.spaceOrder();
+
+ switch (_data.elementType()) {
case TYPE_PNT :
_numLagCoeff = 1;
_dimSimplex = 0;
@@ -550,64 +577,94 @@ void bezierBasis::_construct(int parentType, int p)
subPoints.push_back(gmshGeneratePointsLine(0));
break;
case TYPE_LIN : {
- _numLagCoeff = _order == 0 ? 1 : 2;
+ _numLagCoeff = order == 0 ? 1 : 2;
_dimSimplex = 0;
- _exponents = gmshGenerateMonomialsLine(_order);
- subPoints = generateSubPointsLine(_order);
+ _exponents = gmshGenerateMonomialsLine(order);
+ subPoints = generateSubPointsLine(order);
break;
}
case TYPE_TRI : {
- _numLagCoeff = _order == 0 ? 1 : 3;
+ _numLagCoeff = order == 0 ? 1 : 3;
_dimSimplex = 2;
- _exponents = gmshGenerateMonomialsTriangle(_order);
- subPoints = generateSubPointsTriangle(_order);
+ _exponents = gmshGenerateMonomialsTriangle(order);
+ subPoints = generateSubPointsTriangle(order);
break;
}
case TYPE_QUA : {
- _numLagCoeff = _order == 0 ? 1 : 4;
+ _numLagCoeff = order == 0 ? 1 : 4;
_dimSimplex = 0;
- _exponents = gmshGenerateMonomialsQuadrangle(_order);
- subPoints = generateSubPointsQuad(_order);
+ _exponents = gmshGenerateMonomialsQuadrangle(order);
+ subPoints = generateSubPointsQuad(order);
break;
}
case TYPE_TET : {
- _numLagCoeff = _order == 0 ? 1 : 4;
+ _numLagCoeff = order == 0 ? 1 : 4;
_dimSimplex = 3;
- _exponents = gmshGenerateMonomialsTetrahedron(_order);
- subPoints = generateSubPointsTetrahedron(_order);
+ _exponents = gmshGenerateMonomialsTetrahedron(order);
+ subPoints = generateSubPointsTetrahedron(order);
break;
}
case TYPE_PRI : {
- _numLagCoeff = _order == 0 ? 1 : 6;
+ _numLagCoeff = order == 0 ? 1 : 6;
_dimSimplex = 2;
- _exponents = gmshGenerateMonomialsPrism(_order);
- subPoints = generateSubPointsPrism(_order);
+ _exponents = gmshGenerateMonomialsPrism(order);
+ subPoints = generateSubPointsPrism(order);
break;
}
case TYPE_HEX : {
- _numLagCoeff = _order == 0 ? 1 : 8;
+ _numLagCoeff = order == 0 ? 1 : 8;
_dimSimplex = 0;
- _exponents = gmshGenerateMonomialsHexahedron(_order);
- subPoints = generateSubPointsHex(_order);
+ _exponents = gmshGenerateMonomialsHexahedron(order);
+ subPoints = generateSubPointsHex(order);
break;
}
default : {
Msg::Error("Unknown function space of parentType %d : "
- "reverting to TET_1", parentType);
- _order = 0;
+ "reverting to TET_1", _data.elementType());
_numLagCoeff = 4;
_dimSimplex = 3;
- _exponents = gmshGenerateMonomialsTetrahedron(_order);
- subPoints = generateSubPointsTetrahedron(_order);
+ _exponents = gmshGenerateMonomialsTetrahedron(0);
+ subPoints = generateSubPointsTetrahedron(0);
break;
}
}
_numDivisions = static_cast<int>(subPoints.size());
fullMatrix<double> bezierPoints = _exponents;
- bezierPoints.scale(1./_order);
+ bezierPoints.scale(1./order);
+
+ matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, order, _dimSimplex);
+ matrixBez2Lag.invert(matrixLag2Bez);
+ subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, order, _dimSimplex);
+}
+
+void bezierBasis::_constructPyr()
+{
+ if (_data.elementType() != TYPE_PYR) {
+ Msg::Fatal("This bezierBasis constructor is for pyramids !");
+ }
+
+ const bool pyr = _data.isPyramidalSpace();
+ const int nij = _data.nij(), nk = _data.nk();
- matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, _order, _dimSimplex);
+ _numLagCoeff = nk == 0 ? 4 : 8;
+ _dimSimplex = 0;
+ gmshGenerateMonomials(_data, _exponents);
+
+ fullMatrix<double> bezierPoints;
+ generateBezierPoints(bezierPoints);
+ matrixBez2Lag = generateBez2LagMatrixPyramid(_exponents, bezierPoints,
+ pyr, nij, nk);
matrixBez2Lag.invert(matrixLag2Bez);
- subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, _order, _dimSimplex);
+ if (pyr) {
+ _numDivisions = 0;
+ }
+ else {
+ std::vector< fullMatrix<double> > subPoints;
+ subPoints = generateSubPointsPyr(nij, nk);
+ _numDivisions = static_cast<int>(subPoints.size());
+ subDivisor = generateSubDivisorPyramid(_exponents, subPoints, matrixLag2Bez,
+ pyr, nij, nk);
+ }
+ return;
}
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
index 1b0af7216d..9ca32760d9 100644
--- a/Numeric/bezierBasis.h
+++ b/Numeric/bezierBasis.h
@@ -10,6 +10,7 @@
#include <vector>
#include "fullMatrix.h"
#include "ElementType.h"
+#include "FuncSpaceData.h"
class MElement;
@@ -17,8 +18,8 @@ class bezierBasis {
private :
// the 'numLagCoeff' first exponents are related to 'Lagrangian' values
int _numLagCoeff;
- int _type, _order;
int _numDivisions, _dimSimplex;
+ const FuncSpaceData _data;
friend class MetricBasis;
fullMatrix<double> _exponents;
@@ -28,22 +29,34 @@ class bezierBasis {
fullMatrix<double> matrixBez2Lag;
fullMatrix<double> subDivisor;
+ void printTag() const {Msg::Info("tagBezier is %d", _data.elementTag());}
+ void printFuncSpace() const {_data.print();}
+
// Constructors
- inline bezierBasis(int tag) {
- _construct(ElementType::ParentTypeFromTag(tag), ElementType::OrderFromTag(tag));
- }
- inline bezierBasis(int parendtType, int order) {
- _construct(parendtType, order);
+ inline bezierBasis(FuncSpaceData data) : _data(data) {
+ if (_data.elementType() == TYPE_PYR)
+ _constructPyr();
+ else
+ _construct();
}
// get methods
inline int getDim() const {return _exponents.size2();}
- inline int getOrder() const {return _order;}
+ inline int getOrder() const {return _data.spaceOrder();}
inline int getDimSimplex() const {return _dimSimplex;}
inline int getNumLagCoeff() const {return _numLagCoeff;}
inline int getNumDivision() const {return _numDivisions;}
inline int getNumSubNodes() const {return subDivisor.size1();}
+ // generate Bezier points
+ void generateBezierPoints(fullMatrix<double>&) const;
+
+ // Subdivide Bezier coefficients
+ void subdivideBezCoeff(const fullMatrix<double> &coeff,
+ fullMatrix<double> &subCoeff) const;
+ void subdivideBezCoeff(const fullVector<double> &coeff,
+ fullVector<double> &subCoeff) const;
+
// Interpolation of n functions on N points :
// coeffs(numCoeff, n) and uvw(N, dim)
// => result(N, n)
@@ -51,9 +64,21 @@ class bezierBasis {
const fullMatrix<double> &uvw,
fullMatrix<double> &result,
bool bezCoord = false) const;
+ void interpolate(const fullVector<double> &coeffs,
+ const fullMatrix<double> &uvw,
+ fullVector<double> &result,
+ bool bezCoord = false) const {
+ int size = uvw.size1();
+ result.resize(size);
+ fullMatrix<double> c(const_cast<double*>(coeffs.getDataPtr()), size, 1);
+ fullMatrix<double> r(const_cast<double*>(result.getDataPtr()), size, 1);
+ interpolate(c, uvw, r, bezCoord);
+ }
private :
- void _construct(int parendtType, int order);
+ void _construct();
+ void _constructPyr();
+ void _FEpoints2BezPoints(fullMatrix<double>&) const;
};
#endif
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index e844ee725b..3c01eaf56f 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -6,10 +6,10 @@
#ifndef _FULL_MATRIX_H_
#define _FULL_MATRIX_H_
-#include <math.h>
-#include <stdio.h>
#include "GmshConfig.h"
#include "GmshMessage.h"
+#include <cmath>
+#include <cstdio>
template <class scalar> class fullMatrix;
@@ -150,6 +150,12 @@ class fullVector
return *this;
}
+ void copy(const fullVector<scalar> &v, int i0, int ni, int desti0)
+ {
+ for(int i = i0, desti = desti0; i < i0 + ni; i++, desti++)
+ (*this)(desti) = v(i);
+ }
+
// set
/**
@param r A vector index between 0 and size() - 1;
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index d4fbbfe384..a443c6187b 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -647,3 +647,20 @@ nodalBasis::nodalBasis(int tag)
}
}
+
+void nodalBasis::getReferenceNodesForBezier(fullMatrix<double> &nodes) const {
+ if (parentType != TYPE_PYR && !serendip) {
+ nodes = points;
+ }
+ else {
+ const bool serendipSpace = false;
+ if (parentType != TYPE_PYR) {
+ FuncSpaceData data(true, type, order, &serendipSpace);
+ gmshGeneratePoints(data, nodes);
+ }
+ else {
+ FuncSpaceData data(true, type, false, order, order, &serendipSpace);
+ gmshGeneratePoints(data, nodes);
+ }
+ }
+}
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index 25b855610b..93ab89b38a 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -15,19 +15,23 @@ class nodalBasis {
bool serendip;
fullMatrix<double> points;
- nodalBasis() {};
+ nodalBasis() {}
nodalBasis(int tag);
virtual ~nodalBasis() {}
virtual int getNumShapeFunctions() const = 0;
+ void getReferenceNodes(fullMatrix<double> &nodes) const {
+ nodes = points;
+ }
+ void getReferenceNodesForBezier(fullMatrix<double> &nodes) const;
// Basis functions & gradients evaluation
virtual void f(double u, double v, double w, double *sf) const = 0;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const = 0;
virtual void df(double u, double v, double w, double grads[][3]) const = 0;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const = 0;
- virtual void ddf(double u, double v, double w, double grads[][3][3]) const {Msg::Fatal("Not implemented");};
- virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");};
+ virtual void ddf(double u, double v, double w, double grads[][3][3]) const {Msg::Fatal("Not implemented");}
+ virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");}
// closures is the list of the nodes of each face, in the order of
// the polynomialBasis of the face; fullClosures is mapping of the
@@ -50,27 +54,21 @@ class nodalBasis {
virtual int getClosureType(int id) const { return closures[id].type; }
virtual const std::vector<int> &getClosure(int id) const { return closures[id]; }
virtual const std::vector<int> &getFullClosure(int id) const { return fullClosures[id]; }
- inline int getClosureId(int iFace, int iSign=1, int iRot=0) const;
+ inline int getClosureId(int iFace, int iSign = 1, int iRot = 0) const;
inline void breakClosureId(int i, int &iFace, int &iSign, int &iRot) const;
};
-
-
inline int nodalBasis::getClosureId(int iFace, int iSign, int iRot) const
{
- return iFace + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot;
+ return iFace + numFaces * (iSign == 1 ? 0 : 1) + 2 * numFaces * iRot;
}
-
-
inline void nodalBasis::breakClosureId(int i, int &iFace, int &iSign, int &iRot) const
{
iFace = i % numFaces;
- i = (i - iFace)/numFaces;
+ i = (i - iFace) / numFaces;
iSign = i % 2;
iRot = (i - iSign) / 2;
}
-
-
#endif
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index 69ac1fed18..a0e76f5180 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -4,6 +4,7 @@
// bugs and problems to the public mailing list <gmsh@geuz.org>.
#include "pointsGenerators.h"
+#include "GmshDefines.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
@@ -14,6 +15,48 @@
// Points Generators
+void gmshGeneratePoints(FuncSpaceData data, fullMatrix<double> &points)
+{
+ switch (data.elementType()) {
+ case TYPE_PNT :
+ points = gmshGeneratePointsLine(0);
+ return;
+ case TYPE_LIN :
+ points = gmshGeneratePointsLine(data.spaceOrder());
+ return;
+ case TYPE_TRI :
+ points = gmshGeneratePointsTriangle(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_QUA :
+ points = gmshGeneratePointsQuadrangle(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_TET :
+ points = gmshGeneratePointsTetrahedron(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_PRI :
+ points = gmshGeneratePointsPrism(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_HEX :
+ points = gmshGeneratePointsHexahedron(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_PYR :
+ points = gmshGeneratePointsPyramidGeneral(data.isPyramidalSpace(),
+ data.nij(),
+ data.nk(),
+ data.spaceIsSerendipity());
+ return;
+ default :
+ Msg::Error("Unknown element type %d (tag %d) for points generation",
+ data.elementType(), data.elementTag());
+ return;
+ }
+}
+
fullMatrix<double> gmshGeneratePointsLine(int order)
{
fullMatrix<double> points = gmshGenerateMonomialsLine(order);
@@ -104,8 +147,8 @@ fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, b
if (points.size1() == 1) return points;
+ const int div = pyr ? nk+nij : std::max(nij, nk);
for (int i = 0; i < points.size1(); ++i) {
- int div = pyr ? nk+nij : std::max(nij, nk);
points(i, 2) = (nk - points(i, 2)) / div;
const double scale = 1. - points(i, 2);
points(i, 0) = scale * (-1 + points(i, 0) * 2. / div);
@@ -116,6 +159,48 @@ fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, b
// Monomials Generators
+void gmshGenerateMonomials(FuncSpaceData data, fullMatrix<double> &monomials)
+{
+ switch (data.elementType()) {
+ case TYPE_PNT :
+ monomials = gmshGenerateMonomialsLine(0);
+ return;
+ case TYPE_LIN :
+ monomials = gmshGenerateMonomialsLine(data.spaceOrder());
+ return;
+ case TYPE_TRI :
+ monomials = gmshGenerateMonomialsTriangle(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_QUA :
+ monomials = gmshGenerateMonomialsQuadrangle(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_TET :
+ monomials = gmshGenerateMonomialsTetrahedron(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_PRI :
+ monomials = gmshGenerateMonomialsPrism(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_HEX :
+ monomials = gmshGenerateMonomialsHexahedron(data.spaceOrder(),
+ data.spaceIsSerendipity());
+ return;
+ case TYPE_PYR :
+ monomials = gmshGenerateMonomialsPyramidGeneral(data.isPyramidalSpace(),
+ data.nij(),
+ data.nk(),
+ data.spaceIsSerendipity());
+ return;
+ default :
+ Msg::Error("Unknown element type %d (tag %d) for monomials generation",
+ data.elementType(), data.elementTag());
+ return;
+ }
+}
+
fullMatrix<double> gmshGenerateMonomialsLine(int order, bool serendip)
{
fullMatrix<double> monomials(order + 1, 1);
@@ -786,7 +871,7 @@ fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(
bool pyr, int nij, int nk, bool forSerendipPoints)
{
if (nij < 0 || nk < 0) {
- Msg::Fatal("Wrong arguments for pyramid's monomials generation !");
+ Msg::Fatal("Wrong arguments for pyramid's monomials generation ! (%d & %d)", nij, nk);
}
if (!pyr && nk > 0 && nij == 0) {
Msg::Error("Wrong argument association for pyramid's monomials generation ! Setting nij to 1");
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index 2d581f2f9d..e0a4fc905e 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -7,7 +7,7 @@
#define POINTSGENERATORS_H
#include "fullMatrix.h"
-
+#include "FuncSpaceData.h"
/*
* Functions to generate point distributions on
* the references elements, for all orders.
@@ -18,6 +18,8 @@
// Points
+void gmshGeneratePoints(FuncSpaceData, fullMatrix<double> &);
+
fullMatrix<double> gmshGeneratePointsLine(int order);
fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip = false);
@@ -33,6 +35,8 @@ fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, b
// Monomial exponents
+void gmshGenerateMonomials(FuncSpaceData, fullMatrix<double> &);
+
fullMatrix<double> gmshGenerateMonomialsLine(int order, bool serendip = false);
fullMatrix<double> gmshGenerateMonomialsTriangle(int order, bool serendip = false);
--
GitLab