Commit 1ca5bfe0 authored by Christophe Geuzaine's avatar Christophe Geuzaine

Merge branch 'perf' into 'master'

Minor cleanup

See merge request !146
parents 5fcd9441 8ae6ac5f
Pipeline #1867 passed with stage
in 62 minutes and 17 seconds
......@@ -146,10 +146,10 @@ private:
public:
MEdgeN() {}
MEdgeN(const std::vector<MVertex*> &v);
inline int getNumVertices() const { return (int)_v.size(); }
inline MVertex *getVertex(int i) const { return _v[i]; }
inline const std::vector<MVertex*> &getVertices() const { return _v; }
inline int getPolynomialOrder() const {return getNumVertices() - 1;}
int getNumVertices() const { return (int)_v.size(); }
MVertex *getVertex(int i) const { return _v[i]; }
const std::vector<MVertex*> &getVertices() const { return _v; }
int getPolynomialOrder() const {return getNumVertices() - 1;}
MEdge getEdge() const;
......
......@@ -15,17 +15,20 @@
void MElementBB(void *a, double *min, double *max)
{
MElement *e = (MElement*)a;
if (e->getPolynomialOrder() == 1) {
MElement *e = static_cast<MElement *>(a);
if(e->getPolynomialOrder() == 1) {
MVertex *v = e->getVertex(0);
min[0] = max[0] = v->x();
min[1] = max[1] = v->y();
min[2] = max[2] = v->z();
for(int i = 1; i < e->getNumVertices(); i++){
for(int i = 1; i < e->getNumVertices(); i++) {
v = e->getVertex(i);
min[0] = std::min(min[0], v->x()); max[0] = std::max(max[0], v->x());
min[1] = std::min(min[1], v->y()); max[1] = std::max(max[1], v->y());
min[2] = std::min(min[2], v->z()); max[2] = std::max(max[2], v->z());
min[0] = std::min(min[0], v->x());
max[0] = std::max(max[0], v->x());
min[1] = std::min(min[1], v->y());
max[1] = std::max(max[1], v->y());
min[2] = std::min(min[2], v->z());
max[2] = std::max(max[2], v->z());
}
}
else {
......@@ -33,23 +36,24 @@ void MElementBB(void *a, double *min, double *max)
e->getNodesCoord(nodesXYZ);
fullMatrix<double> bezNodes(e->getNumVertices(), 3);
const bezierBasis *bez = BasisFactory::getBezierBasis(FuncSpaceData(e));
bez->lag2Bez(nodesXYZ, bezNodes);
min[0] = max[0] = bezNodes(0, 0);
min[1] = max[1] = bezNodes(0, 1);
min[2] = max[2] = bezNodes(0, 2);
for(int i = 1; i < e->getNumVertices(); i++){
min[0] = std::min(min[0], bezNodes(0, 0));
max[0] = std::max(max[0], bezNodes(0, 0));
min[1] = std::min(min[1], bezNodes(0, 1));
max[1] = std::max(max[1], bezNodes(0, 1));
min[2] = std::min(min[2], bezNodes(0, 2));
max[2] = std::max(max[2], bezNodes(0, 2));
for(int i = 1; i < e->getNumVertices(); i++) {
min[0] = std::min(min[0], bezNodes(i, 0));
max[0] = std::max(max[0], bezNodes(i, 0));
min[1] = std::min(min[1], bezNodes(i, 1));
max[1] = std::max(max[1], bezNodes(i, 1));
min[2] = std::min(min[2], bezNodes(i, 2));
max[2] = std::max(max[2], bezNodes(i, 2));
}
}
// make bounding boxes larger up to (absolute) geometrical tolerance
double eps = CTX::instance()->geom.tolerance;
for(int i = 0; i < 3; i++){
double const eps = CTX::instance()->geom.tolerance;
for(int i = 0; i < 3; i++) {
min[i] -= eps;
max[i] += eps;
}
......
......@@ -14,10 +14,9 @@
#include "qualityMeasures.h"
#endif
#include <cmath>
#include <cstring>
#define SQU(a) ((a)*(a))
void MTriangle::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
SVector3 *n)
{
......@@ -72,20 +71,20 @@ double MTriangle::getInnerRadius()
dist[i] = e.getVertex(0)->distance(e.getVertex(1));
k += 0.5 * dist[i];
}
double area = sqrt(k*(k-dist[0])*(k-dist[1])*(k-dist[2]));
double const area = std::sqrt(k*(k-dist[0])*(k-dist[1])*(k-dist[2]));
return area/k;
}
double MTriangle::getOuterRadius()
{
// radius of circle circumscribing a triangle
double dist[3], k = 0.;
double dist[3], k = 0.0;
for (int i = 0; i < 3; i++){
MEdge e = getEdge(i);
dist[i] = e.getVertex(0)->distance(e.getVertex(1));
k += 0.5 * dist[i];
}
double area = sqrt(k*(k-dist[0])*(k-dist[1])*(k-dist[2]));
double const area = std::sqrt(k*(k-dist[0])*(k-dist[1])*(k-dist[2]));
return dist[0]*dist[1]*dist[2]/(4*area);
}
......@@ -132,20 +131,21 @@ void MTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
R[1] = (xyz[1] - p0.y());
sys2x2(M, R, uvw);
return;*/
const double O[3] = {_v[0]->x(), _v[0]->y(), _v[0]->z()};
const double d[3] = {xyz[0] - O[0], xyz[1] - O[1], xyz[2] - O[2]};
const double d1[3] = {_v[1]->x() - O[0], _v[1]->y() - O[1], _v[1]->z() - O[2]};
const double d2[3] = {_v[2]->x() - O[0], _v[2]->y() - O[1], _v[2]->z() - O[2]};
const double Jxy = d1[0] * d2[1] - d1[1] * d2[0];
const double Jxz = d1[0] * d2[2] - d1[2] * d2[0];
const double Jyz = d1[1] * d2[2] - d1[2] * d2[1];
if ((fabs(Jxy) > fabs(Jxz)) && (fabs(Jxy) > fabs(Jyz))){
if ((std::abs(Jxy) > std::abs(Jxz)) && (std::abs(Jxy) > std::abs(Jyz))){
uvw[0] = (d[0] * d2[1] - d[1] * d2[0]) / Jxy;
uvw[1] = (d[1] * d1[0] - d[0] * d1[1]) / Jxy;
}
else if (fabs(Jxz) > fabs(Jyz)){
else if (std::abs(Jxz) > std::abs(Jyz)){
uvw[0] = (d[0] * d2[2] - d[2] * d2[0]) / Jxz;
uvw[1] = (d[2] * d1[0] - d[0] * d1[2]) / Jxz;
}
......@@ -153,7 +153,7 @@ void MTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
uvw[0] = (d[1] * d2[2] - d[2] * d2[1]) / Jyz;
uvw[1] = (d[2] * d1[1] - d[1] * d1[2]) / Jyz;
}
uvw[2] = 0.;
uvw[2] = 0.0;
}
int MTriangle::numCommonNodesInDualGraph(const MElement *const other) const
......@@ -237,11 +237,11 @@ bool MTriangle::getFaceInfo(const MFace & face, int &ithFace, int &sign,
int MTriangle6::getNumFacesRep(bool curved)
{
return curved ? SQU(CTX::instance()->mesh.numSubEdges) : 1;
return curved ? std::pow(CTX::instance()->mesh.numSubEdges, 2) : 1;
}
int MTriangleN::getNumFacesRep(bool curved)
{
return curved ? SQU(CTX::instance()->mesh.numSubEdges) : 1;
return curved ? std::pow(CTX::instance()->mesh.numSubEdges, 2) : 1;
}
static void _myGetFaceRep(MTriangle *t, int num, double *x, double *y, double *z,
......
......@@ -289,12 +289,12 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
bool Extend1dMeshIn2dSurfaces()
{
return CTX::instance()->mesh.lcExtendFromBoundary ? true : false;
return CTX::instance()->mesh.lcExtendFromBoundary;
}
bool Extend2dMeshIn3dVolumes()
{
return CTX::instance()->mesh.lcExtendFromBoundary ? true : false;
return CTX::instance()->mesh.lcExtendFromBoundary;
}
SMetric3 max_edge_curvature_metric(const GVertex *gv)
......
......@@ -35,7 +35,7 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
extern double lengthInfniteNorm(const double p[2], const double q[2],
const double quadAngle);
inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f)
double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u),
0.5 * (p1->v + p2->v)));
......@@ -46,28 +46,28 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f)
const double dx1 = p1->X - GP.x();
const double dy1 = p1->Y - GP.y();
const double dz1 = p1->Z - GP.z();
const double l1 = sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
const double l1 = std::sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
const double dx2 = p2->X - GP.x();
const double dy2 = p2->Y - GP.y();
const double dz2 = p2->Z - GP.z();
const double l2 = sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
const double l2 = std::sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
return l1 + l2;
}
inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f)
double computeEdgeLinearLength(BDS_Edge *e, GFace *f)
{
// FIXME !!!
if (f->geomType() == GEntity::Plane)
return e->length();
else
return computeEdgeLinearLength(e->p1, e->p2, f);
return f->geomType() == GEntity::Plane ?
e->length() :
computeEdgeLinearLength(e->p1, e->p2, f);
}
double NewGetLc(BDS_Point *p)
double NewGetLc(BDS_Point *point)
{
return Extend1dMeshIn2dSurfaces() ?
std::min(p->lc(), p->lcBGM()) : p->lcBGM();
return Extend1dMeshIn2dSurfaces() ? std::min(point->lc(), point->lcBGM()) :
point->lcBGM();
}
static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f)
......@@ -86,7 +86,7 @@ static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f)
if(CTX::instance()->mesh.lcFromCurvature){
double l3 = l;
double lcmin = std::min(std::min(l1, l2), l3);
double const lcmin = std::min(std::min(l1, l2), l3);
l1 = std::min(lcmin*1.2,l1);
l2 = std::min(lcmin*1.2,l2);
l3 = std::min(lcmin*1.2,l3);
......@@ -95,12 +95,10 @@ static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f)
return l;
}
double NewGetLc(BDS_Edge *e, GFace *f)
double NewGetLc(BDS_Edge *const edge, GFace *const face)
{
double linearLength = computeEdgeLinearLength(e, f);
double l = correctLC_ (e->p1,e->p2,f);
// printf("BDS %d %d %g %g correct lc =%g lreal=%g \n",e->p1->iD, e->p2->iD,e->p1->lc(),e->p2->lc(),l,linearLength);
return linearLength / l;
return computeEdgeLinearLength(edge, face) /
correctLC_(edge->p1, edge->p2, face);
}
double NewGetLc(BDS_Point *p1, BDS_Point *p2, GFace *f)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment