From a6d1f76e2a017b9131cce9a7596c18ed435abfb1 Mon Sep 17 00:00:00 2001
From: Gaetan Bricteux <gaetan.bricteux@uclouvain.be>
Date: Thu, 6 May 2010 09:20:29 +0000
Subject: [PATCH] fix MLine::getIntegrationPoints
---
Geo/MElementCut.cpp | 44 +++++++++----------
Geo/MLine.cpp | 14 +-----
contrib/DiscreteIntegration/Integration3D.cpp | 2 +-
3 files changed, 23 insertions(+), 37 deletions(-)
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index f16e3bd60a..eb4267ecf0 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -267,14 +267,10 @@ bool MLineChild::isInside(double u, double v, double w)
void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
-
-
*npts = 0;
double jac[3][3];
if(_intpt) delete [] _intpt;
-
- int nbP = pOrder / 2 + 1; // MLine::getIntegrationPoints()
- _intpt = new IntPt[nbP];
+ _intpt = new IntPt[getNGQLPts(pOrder)];
int nptsi;
IntPt *ptsi;
@@ -291,30 +287,30 @@ void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
// -----------mich mach---------------------//
- MVertex *vo = _orig->getVertex(0);
- MVertex *vf = _orig->getVertex(1);
+ MVertex *vo = _orig->getVertex(0);
+ MVertex *vf = _orig->getVertex(1);
+
+ SPoint3 P(vo->x(), vo->y(), vo->z());
+ SPoint3 Q(vf->x(), vf->y(), vf->z());
+
+ SPoint3 R;
+
+ R = P + Q;
+ R/=2;
- SPoint3 P(vo->x(), vo->y(), vo->z());
- SPoint3 Q(vf->x(), vf->y(), vf->z());
+ vo = getVertex(0);
+ vf = getVertex(1);
- SPoint3 R;
-
- R = P + Q;
- R/=2;
+ SPoint3 A(vo->x(), vo->y(), vo->z());
+ SPoint3 B(vf->x(), vf->y(), vf->z());
- vo = getVertex(0);
- vf = getVertex(1);
+ double lengthDemi = R.distance(Q);
- SPoint3 A(vo->x(), vo->y(), vo->z());
- SPoint3 B(vf->x(), vf->y(), vf->z());
+ if (P.distance(A) < Q.distance(A)) {v_uvw[0][0] = - R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
+ else {v_uvw[0][0] = R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
- double lengthDemi = R.distance(Q);
-
- if (P.distance(A) < Q.distance(A)) {v_uvw[0][0] = - R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
- else {v_uvw[0][0] = R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
-
- if (P.distance(B) < Q.distance(B)) {v_uvw[1][0] = - R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
- else {v_uvw[1][0] = R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
+ if (P.distance(B) < Q.distance(B)) {v_uvw[1][0] = - R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
+ else {v_uvw[1][0] = R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
// ---------------------------------------//
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index b6cdc2fc7f..eb79c54a6c 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -31,18 +31,8 @@ const polynomialBasis* MLine::getFunctionSpace(int o) const
void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
- static IntPt GQL[100];
- double *t, *w;
- int nbP = pOrder / 2 + 1;
- gmshGaussLegendre1D(nbP, &t, &w);
- for (int i = 0; i < nbP; i++){
- GQL[i].pt[0] = t[i];
- GQL[i].pt[1] = 0;
- GQL[i].pt[2] = 0;
- GQL[i].weight = w[i];
- }
- *npts = nbP;
- *pts = GQL;
+ *npts = getNGQLPts(pOrder);
+ *pts = getGQLPts(pOrder);
}
double MLine::getInnerRadius()
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 4dd2e58936..48b7bbf9b1 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -3,7 +3,7 @@
#include "Integration3D.h"
#include "recurCut.h"
-#include "../../Numeric/Gauss.h"
+#include "Gauss.h"
#define ZERO_LS_TOL 1.e-7
#define EQUALITY_TOL 1.e-15
--
GitLab