From e917b465c3b9e69ce7b58bd74691357262d6aea0 Mon Sep 17 00:00:00 2001
From: SYLVAND Guillaume <guillaume.sylvand@airbus.com>
Date: Mon, 10 Feb 2025 18:25:02 +0100
Subject: [PATCH] sync with test_FEMBEM@airbus + cosmetics
- different way to create pipe mesh to have a correct (=conformal) mesh
Note: it changes (lowers) the number of non-zeros in the FEM matrices
---
include/main.h | 7 +++-
include/mpf_simd.h | 4 +++
src/cylinder.c | 75 +++++++++++++++++++++-----------------
src/hchameleon.c | 2 ++
src/help.c | 3 +-
src/kernel.c | 6 ++--
src/main.c | 12 ++++++-
src/pipe.c | 89 ++++++++++++++++++++++++++++------------------
src/prepareTEST.c | 2 +-
src/rhs.c | 5 +--
src/testHMAT.c | 2 +-
11 files changed, 131 insertions(+), 76 deletions(-)
create mode 100644 include/mpf_simd.h
diff --git a/include/main.h b/include/main.h
index aab2feb..5716f25 100644
--- a/include/main.h
+++ b/include/main.h
@@ -12,6 +12,7 @@
#endif
/* ================================================================================= */
#include <math.h>
+#include <complex.h>
#include "mpi.h"
#ifdef HAVE_HMAT
#include "hmat/hmat.h"
@@ -127,6 +128,11 @@ extern int writeMesh;
/*! \brief Write a UNV file of the mesh. Enabled with option --write-mesh-unv. */
extern int writeMeshUNV;
+/*! \brief Create a conformal volume mesh by using alternatively 2 ways of splitting hexa into 5 tetra.
+
+ Disabled with --no-conformal-mesh, enabled with --conformal-mesh. Default is enabled. */
+extern int conformalMesh;
+
/*! \brief List of algorithms that we can test. */
enum algo {
_undefined,
@@ -246,7 +252,6 @@ int computeRelativeError(void *sol, void *ref, double *eps) ;
int computeVecNorm(void *ref, double *norm);
int sommePara(double *buf) ;
int computeCoordCylinder(int i, double *coord) ;
-int computeCoordCylinderDetail(int i, double *coord) ;
int initCylinder(int *argc, char ***argv) ;
int initCylinderDetail(int *argc, char ***argv) ;
int getMeshStep(double *m) ;
diff --git a/include/mpf_simd.h b/include/mpf_simd.h
new file mode 100644
index 0000000..7889386
--- /dev/null
+++ b/include/mpf_simd.h
@@ -0,0 +1,4 @@
+// Replacement for MPF macros & routines for SIMD optimisations
+#define MPF_SIMD_DISP(x) x
+#define mpf_simd_cexp cexp
+#define mpf_simd_sincos(theta,s,c) {*s=sin(theta);*(c)=cos(theta);}
diff --git a/src/cylinder.c b/src/cylinder.c
index d50087c..ccc9e67 100644
--- a/src/cylinder.c
+++ b/src/cylinder.c
@@ -1,6 +1,5 @@
#include "main.h"
-
-/* ============================================================================================================ */
+#include <mpf_simd.h>
/*! \brief Radius of the cylinder */
double radius = 2.0 ;
@@ -49,40 +48,25 @@ inline static int computeCoordCylinderMain(int i, double *coord) {
/* zStep and angleStep are defined by :
angleStep*radius = zStep*nbPtsLoop (= meshStep)
angleStep*nbPtsLoop = 2.pi */
- zStep=height/(double)nbPtsMain ;
+ zStep=height/(double)(nbPtsMain-1) ;
angleStep=sqrt(zStep*2.*M_PI/radius) ;
meshStep=zStep*2.*M_PI/angleStep ;
nbPtsLoop=floor(2.*M_PI/angleStep) ; // Used to write VTK file
+ // Once I have 'nbPtsLoop' as an integer, I recompute meshStep and angleStep
+ angleStep = 2.*M_PI/(double)nbPtsLoop;
+ meshStep = angleStep*radius;
}
theta = (double)i * angleStep ;
- coord[0] = radius * sin(theta) ;
- coord[1] = radius * cos(theta) ;
+ mpf_simd_sincos(theta, coord, coord+1);
+ coord[0] *= radius;
+ coord[1] *= radius;
coord[2] = (double)i*zStep ;
return 0 ;
}
-/*! \brief Computes coordinates of point number 'i' on the object (main cylinder or detail).
-
- \a i should be in [0, nbPts[.
- \param i the number of the point in the detail (input)
- \param coord the coordinates x y z of this point (output)
- \return 0 for success
-*/
-int computeCoordCylinder(int i, double *coord) {
- ASSERTQ(i<nbPtsBEM && i>=0);
- if (testedModel == _fembem)
- return computeCoordPipe(i, coord);
- if (useDetail && i >= nbPtsMain) {
- i = i - nbPtsMain;
- return computeCoordCylinderDetail(i, coord);
- } else {
- return computeCoordCylinderMain(i, coord);
- }
-}
-
/*! \brief Computes coordinates of point number 'i' of the detail placed on the main cylinder.
The points are computed on a cylinder defined by 'heightDetail' and 'radiusDetail' (values in number of wavelength).
@@ -93,39 +77,63 @@ int computeCoordCylinder(int i, double *coord) {
\param coord the coordinates x y z of this point (output)
\return 0 for success
*/
-int computeCoordCylinderDetail(int i, double *coord) {
+inline static int computeCoordCylinderDetail(int i, double *coord) {
ASSERTQ(testedModel == _bem);
static double angleStep = 0. ;
static double zStep = 0. ;
+ static double heightD = 0.;
+ static double radiusD = 0.;
double theta ;
- double heightD = lambda * heightDetail ;
- double radiusD = lambda * radiusDetail ;
-
if (i<0 || i>=(nbPts - nbPtsMain)) {
SETERRQ(1, "Incorrect unknown number %d, should be in [0, %d[\n", i, nbPts - nbPtsMain) ;
}
if (angleStep == 0.) {
+ heightD = lambda * heightDetail ;
+ radiusD = lambda * radiusDetail ;
/* zStep and angleStep are defined by :
angleStep*radius = zStep*nbPtsLoopDetail (= meshStep)
angleStep*nbPtsLoopDetail = 2.pi */
- zStep = heightD/(double)(nbPts - nbPtsMain) ;
+ zStep = heightD/(double)(nbPts - nbPtsMain-1) ;
angleStep = sqrt(zStep*2.*M_PI/radiusD) ;
meshStepDetail = zStep*2.*M_PI/angleStep ;
nbPtsLoopDetail=floor(2.*M_PI/angleStep) ; // Used to write VTK file
+ // Once I have 'nbPtsLoopDetail' as an integer, I recompute meshStepDetail and angleStep
+ angleStep = 2.*M_PI/(double)nbPtsLoopDetail;
+ meshStepDetail = angleStep*radius;
}
/* Le detail est oriente selon y, perpendiculaire au cylindre principal, a mi-hauteur */
theta = (double)i * angleStep ;
- coord[0] = radiusD * sin(theta) ;
- coord[1] = zStep * (double)(i+1) + radius ;
- coord[2] = radiusD * cos(theta) + 0.5 * height ;
+ mpf_simd_sincos(theta, coord, coord+2);
+ coord[0] *= radiusD ;
+ coord[1] = zStep * (double)(i+1) + radius ; // 'i+1' instead of 'i' to avoid that the detail touches the main cylinder
+ coord[2] = radiusD * coord[2] + 0.5 * height ;
return 0 ;
}
+/*! \brief Computes coordinates of point number 'i' on the object (main cylinder or detail).
+
+ \a i should be in [0, nbPts[.
+ \param i the number of the point in the detail (input)
+ \param coord the coordinates x y z of this point (output)
+ \return 0 for success
+*/
+MPF_SIMD_DISP(int computeCoordCylinder(int i, double *coord)) {
+ ASSERTQ(i<nbPtsBEM && i>=0);
+ if (testedModel == _fembem)
+ return computeCoordPipe(i, coord);
+ if (useDetail && i >= nbPtsMain) {
+ i = i - nbPtsMain;
+ return computeCoordCylinderDetail(i, coord);
+ } else {
+ return computeCoordCylinderMain(i, coord);
+ }
+}
+
double* createCylinder(void) {
double* result = MpfCalloc((size_t)3 * nbPts, sizeof(double));
int i, ierr;
@@ -325,8 +333,9 @@ int initCylinderDetail(int *argc, char ***argv) {
/* Nombre de points sur le detail = superficie du cylindre (2.pi.R.h.lambda^2) divise par superficie d'une "maille" (step^2) */
nbPts += ceil(2*M_PI*lambda*lambda*heightDetail*radiusDetail/(step*step));
+ nbPtsBEM = nbPts ;
- /* Calcule les coordonnees d'un point, afin de fixer meshStepDetail (qui peut servir a fixer lambda) */
+ /* Calcule les coordonnees d'un point, afin de fixer meshStepDetail */
ierr = computeCoordCylinder(nbPtsMain, &(coord[0]) ) ; CHKERRQ(ierr) ;
return 0 ;
diff --git a/src/hchameleon.c b/src/hchameleon.c
index f812833..04c6c83 100644
--- a/src/hchameleon.c
+++ b/src/hchameleon.c
@@ -1,3 +1,5 @@
+#include "config.h"
+
#if defined(HAVE_MKL_H)
#include <mkl.h>
#else
diff --git a/src/help.c b/src/help.c
index 9ed0ba4..41006ee 100644
--- a/src/help.c
+++ b/src/help.c
@@ -10,6 +10,7 @@ int printHelp() {
" -nbpts: Number of unknowns of the global linear system. Default: 16000 (can be in floating point notation, like 1.2e4).\n"
" -nbrhs: Number of right-hand sides of the global linear system. Default: 1.\n"
" -s/-d/-c/-z: Scalar type of the computation (real computation uses 1/R kernel in BEM, complex computation uses oscillatory exp(ikr)/r kernel). Default: z.\n"
+ " --partial-fembem: Using partial interactions between BEM and FEM unknowns. Default is NO."
" -v: Verbose mode, to display values of data (extremely verbose). Default is OFF.\n"
" -sparserhs: Sparsity of the right hand sides (1 non-zero value every 'n' values). Default is nbpts/80.log10(nbpts).\n"
" --new-rhs/--no-new-rhs: Use (or not) new way to compute RHS (old way was producing low rank RHS). Default is ON.\n"
@@ -24,7 +25,7 @@ int printHelp() {
" -height: Height of the cylinder. Default is 4.\n"
" -with_detail: add a finely meshed cylionder detail. Default is NO.\n"
" -radius_detail: Radius of the cylinder detail. Default is 0.02.\n"
- " -height: Height of the cylinder detail. Default is 0.5.\n"
+ " -height_detail: Height of the cylinder detail. Default is 0.5.\n"
" -ptsperlambda_detail: Points per wavelength on the detail. Default is '10 times finer than on the main cylinder'.\n"
"\n"
" FEM and FEM/BEM tests use a volumic cylinder, defined with:\n"
diff --git a/src/kernel.c b/src/kernel.c
index 5415bd0..ae28e7a 100644
--- a/src/kernel.c
+++ b/src/kernel.c
@@ -1,6 +1,8 @@
#include "main.h"
+#include <mpf_simd.h>
+
extern double meshStep;
-int computeKernelBEM(double *coord1, double *coord2, int self, double complex *kernel) {
+MPF_SIMD_DISP(int computeKernelBEM(double *coord1, double *coord2, int self, double complex *kernel)) {
double r, v[3], k=0. ;
v[0]=coord1[0]-coord2[0] ;
@@ -11,7 +13,7 @@ int computeKernelBEM(double *coord1, double *coord2, int self, double complex *k
if (complexALGO) {
k=2*M_PI/lambda;
- *kernel = cexp(I * k * r) / (4 * M_PI * r);
+ *kernel = mpf_simd_cexp(I * k * r) / (4 * M_PI * r);
} else
*kernel = 1 / (4 * M_PI * r);
diff --git a/src/main.c b/src/main.c
index 0d6efda..6182774 100644
--- a/src/main.c
+++ b/src/main.c
@@ -38,6 +38,7 @@ double ptsPerLambdaDetail = -1.;
double radiusLeaf = 0.;
int writeMesh = 0;
int writeMeshUNV = 0;
+int conformalMesh=1;
enum algo testedAlgo = _undefined;
int use_simple_assembly = 0;
int divider = 2;
@@ -221,6 +222,15 @@ int main(int argc, char **argv) {
writeMeshUNV=1;
Mpf_printf(MPI_COMM_WORLD, "Activate writing of UNV file 'mesh.unv'\n") ;
}
+ /* --- Create (or not) a conformal volume mesh in pipe.c --- */
+ if (MpfArgHasName(&argc, argv, 1, "--conformal-mesh")) {
+ conformalMesh=1;
+ Mpf_printf(MPI_COMM_WORLD, "Create a conformal volume mesh.\n") ;
+ }
+ if (MpfArgHasName(&argc, argv, 1, "--no-conformal-mesh")) {
+ conformalMesh=0;
+ Mpf_printf(MPI_COMM_WORLD, "Do not create a conformal volume mesh (old way)\n") ;
+ }
switch(testedModel){
case _bem: // 2D mesh
@@ -254,7 +264,7 @@ int main(int argc, char **argv) {
if (MpfArgHasName(&argc, argv, 1, "-with_detail")) {
useDetail = 1;
ASSERTQ(testedModel==_bem); // detail is not (yet) compatible with FEM
- Mpf_printf(MPI_COMM_WORLD, "Add an overmeshed detail");
+ Mpf_printf(MPI_COMM_WORLD, "Add an overmeshed detail\n");
}
/* Wavelength */
diff --git a/src/pipe.c b/src/pipe.c
index b96e6f7..0041763 100644
--- a/src/pipe.c
+++ b/src/pipe.c
@@ -50,9 +50,10 @@ int computeCoordPipe(int i, double *coord) {
if (angleStep==0.) {
/* The pipe is defined as follows: its height is 'height', the points are located between 'radius'
- and 'radius/2'. We have n_r, n_z, n_l points in the direction z, radial and orthoradial.
+ and 'radius/2'. We have n_r, n_z, n_l points in the direction radial, z and orthoradial.
total number of points nbPts=n_r.n_z.n_l
2.pi*radius/n_l = heigth/n_z = (radius/2)/n_r = meshStep
+ (rigorously, we should solve: 2.pi*radius/n_l = heigth/(n_z-1) = (radius/2)/(n_r-1) = meshStep)
We compute the number of points as double and round them up to integers.
angleStep*n_l = 2.pi
zStep=heigth/(n_z.n_l)
@@ -61,8 +62,10 @@ int computeCoordPipe(int i, double *coord) {
double dn_r=radius*dn_z/2./height;
double dn_l=4.*M_PI*dn_r;
n_l=(int)ceil(dn_l);
+ if (conformalMesh && n_l%2 == 0) n_l++; // To create a conformal mesh, n_l must be odd
n_r=(int)ceil(dn_r);
n_z=(int)ceil(dn_z);
+ while(n_l*n_r*n_z < nbPts) n_z++;
meshStep=height/(double)n_z;
zStep=meshStep/(double)n_l;
angleStep=2.*M_PI/(double)n_l;
@@ -142,11 +145,15 @@ int freeFemMatrix() {
}
int createFemMatrix() {
- int nbHexa=0, ierr, rank, iVtk;
- int defTetra[5][4]= { {0,1,2,4}, {1,4,5,7}, {1,2,3,7}, {2,4,6,7}, {1,2,4,7}};
+ int ierr, rank, iVtk;
+ // in a conformal mesh, we alternatively have 'A' and 'B' hexa to match diagonal edges between consecutive hexa
+ // in a non-conformal mesh (the old way), we use only 'A' hexa
+ // we can switch between the 2 types with --conformal-mesh and --no-conformal-mesh, conformal is the default
+ int defTetra[2][5][4]= { { {0,1,2,4}, {1,4,5,7}, {1,2,3,7}, {2,4,6,7}, {1,2,4,7}}, // 'A' hexa: face 0246 is split along 24
+ { {0,1,3,5}, {0,2,3,6}, {0,4,5,6}, {3,5,6,7}, {0,3,5,6}} }; // 'B' hexa: face 0246 is split along 06
FILE *fvtk=NULL, *funv=NULL;
int nbTetra=0, nbTria=0;
-
+ fpos_t posNbTetra;
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if (rank==0) {
if (writeMesh) {
@@ -156,14 +163,16 @@ int createFemMatrix() {
fprintf(fvtk, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
"byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n");
fprintf(fvtk, "\t<UnstructuredGrid>\n");
- fprintf(fvtk, "\t\t<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",
- n_r*n_z*n_l, 5*(n_r*n_l*n_z-n_l*n_z-(n_r-1)*(n_l+1)));
+ fprintf(fvtk, "\t\t<Piece NumberOfPoints=\"%d\" NumberOfCells=\"", nbPts);
+ fgetpos(fvtk, &posNbTetra); // we will write the exact number of tetra later, within the empty space left below
+ fprintf(fvtk, " >\n");
+
// Write the vertices
fprintf(fvtk, "\t\t\t<Points>\n");
fprintf(fvtk, "\t\t\t\t<DataArray type=\"Float32\""
" NumberOfComponents=\"3\" format=\"ascii\">\n");
double coord[3];
- for (iVtk = 0; iVtk < n_r*n_z*n_l; iVtk++) {
+ for (iVtk = 0; iVtk < nbPts; iVtk++) {
ierr = computeCoordPipe(iVtk, coord); CHKERRQ(ierr);
fprintf(fvtk, "%g %g %g\n", (float) (coord[0]), (float)(coord[1]), (float)(coord[2]));
}
@@ -191,7 +200,7 @@ int createFemMatrix() {
fprintf(funv, "%s\n",sUNV_SEPARATOR ) ;
fprintf(funv, "%s\n", sNODE_UNV_ID );
double coord[3];
- for(int i = 1; i<=n_r*n_z*n_l; i++ ) {
+ for(int i = 1; i<=nbPts; i++ ) {
ierr = computeCoordPipe(i-1, coord); CHKERRQ(ierr);
fprintf(funv, "%10d%s\n", i, sNODE_UNV_DESCR );
fprintf(funv, "%25.16E%25.16E%25.16E\n", coord[0], coord[1], coord[2]) ;
@@ -219,10 +228,10 @@ int createFemMatrix() {
// the outer face is based on points i_face, i_face+1, i_face+n_l, i_face+n_l+1, they must be in [0, n_z.n_l-1]
if (i_face+n_l+1 > n_z*n_l-1) continue;
- nbHexa++;
+ // In a conformal mesh, we use alternatively the 2 types of hexa defined in defTetra[]
+ int typeHexa = conformalMesh ? (i_face+i_helix)%2 : 0;
- /* Vertices of the current hexaedre, in [0..nbPts[
- */
+ /* Vertices of the current hexaedre, in [0..nbPts[ */
int i_vert[8];
i_vert[0] = i_helix * n_z*n_l + i_face; // actually it is i_hexa
i_vert[1] = i_helix * n_z*n_l + i_face+1;
@@ -240,34 +249,33 @@ int createFemMatrix() {
int i_vert_t[4]; // vertices of the current tetraedron, in [0..nbPts[
int i, j;
for (i=0 ; i<4 ; i++)
- i_vert_t[i] = i_vert[defTetra[i_tetra][i]];
+ i_vert_t[i] = i_vert[defTetra[typeHexa][i_tetra][i]];
// Loop on the unkowns (=vertices) in this tetrahedron
// Each unknown interacts with the other unknown in the same tetrahedron
// Value of the interaction is 0.1 for self, 0.05 otherwise
- for (i=0 ; i<4 ; i++)
- if (i_vert_t[i] >= 0 && i_vert_t[i] < nbPts)
+ if (i_vert_t[0] >= 0 && i_vert_t[0] < nbPts && i_vert_t[1] >= 0 && i_vert_t[1] < nbPts && i_vert_t[2] >= 0 && i_vert_t[2] < nbPts && i_vert_t[3] >= 0 && i_vert_t[3] < nbPts) {
+ nbTetra++;
+ for (i=0 ; i<4 ; i++)
for (j=0 ; j<4 ; j++)
- if (i_vert_t[j] >= 0 && i_vert_t[j] < nbPts)
- addValue(i_vert_t[i], i_vert_t[j], i==j ? 0.1 : 0.05);
-
- if (writeMesh)
- fprintf(fvtk, "%d %d %d %d\n", i_vert_t[0], i_vert_t[1], i_vert_t[2], i_vert_t[3]);
- if (writeMeshUNV) {
- fprintf(funv, "%10d%s\n", ++nbTetra, sELT_TETRA4_DESC ) ;
- fprintf(funv, "%10d%10d%10d%10d\n", i_vert_t[0]+1, i_vert_t[1]+1, i_vert_t[2]+1, i_vert_t[3]+1);
+ addValue(i_vert_t[i], i_vert_t[j], i==j ? 0.1 : 0.05);
+
+ if (writeMesh)
+ fprintf(fvtk, "%d %d %d %d\n", i_vert_t[0], i_vert_t[1], i_vert_t[2], i_vert_t[3]);
+ if (writeMeshUNV) {
+ fprintf(funv, "%10d%s\n", nbTetra, sELT_TETRA4_DESC ) ;
+ fprintf(funv, "%10d%10d%10d%10d\n", i_vert_t[0]+1, i_vert_t[1]+1, i_vert_t[2]+1, i_vert_t[3]+1);
+ }
}
}
}
if (writeMesh) {
- // Check that the number of hexas created matches the one computed at the beginning and written in the vtk header
- ASSERTQ(nbHexa == n_r*n_l*n_z-n_l*n_z-(n_r-1)*(n_l+1));
fprintf(fvtk, "\t\t\t\t</DataArray>\n");
/* Ecriture des offsets */
fprintf(fvtk, "\t\t\t\t<DataArray "
"type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n");
- for (iVtk = 0; iVtk < nbHexa*5; iVtk++) {
+ for (iVtk = 0; iVtk < nbTetra; iVtk++) {
fprintf(fvtk, "%d ", iVtk*4+4);
if (iVtk % 10 == 9)
fprintf(fvtk, "\n");
@@ -277,7 +285,7 @@ int createFemMatrix() {
/* Ecriture des types d'elements */
fprintf(fvtk, "\t\t\t\t<DataArray "
"type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
- for (iVtk = 0; iVtk < nbHexa*5; iVtk++) {
+ for (iVtk = 0; iVtk < nbTetra; iVtk++) {
fprintf(fvtk, "%d ", 10); // For VTK_TETRA
if (iVtk % 10 == 9)
fprintf(fvtk, "\n");
@@ -285,27 +293,40 @@ int createFemMatrix() {
fprintf(fvtk, "\t\t\t\t</DataArray>\n");
fprintf(fvtk, "\t\t\t</Cells>\n");
fprintf(fvtk, "\t\t</Piece>\n\t</UnstructuredGrid>\n</VTKFile>\n");
+
+ fsetpos(fvtk, &posNbTetra); // Write a-posteriori the exact number of hexa
+ fprintf(fvtk, "%d\"", nbTetra) ;
+
fclose(fvtk);
Mpf_printf(MPI_COMM_WORLD, "Done writing 'mesh.vtu'\n") ;
}
if (writeMeshUNV) {
// Add the triangles on the outer helix
+ // Like for tetra, the shape of triangles depends on conformal/non-conformal mesh
+ int defTria[2][2][3]= { { {0,1,2}, {2,1,3}}, // 'A' hexa
+ { {0,1,3}, {2,0,3}} }; // 'B' hexa
for (int i_face=0 ; i_face<n_z*n_l ; i_face++) {
// the outer face is based on points i_face, i_face+1, i_face+n_l, i_face+n_l+1, they must be in [0, n_z.n_l-1]
if (i_face+n_l+1 > n_z*n_l-1) continue;
- /* Vertices of the current outer triangles in [0..n_z*n_l[
- */
+ // In a conformal mesh, we use alternatively the 2 types of hexa defined in defTetra[]
+ int typeHexa = conformalMesh ? i_face%2 : 0;
+
+ /* Vertices of the current outer triangles in [0..n_z*n_l[ */
int i_vert[4];
i_vert[0] = i_face;
i_vert[1] = i_face+1;
i_vert[2] = i_face+n_l;
i_vert[3] = i_face+n_l+1;
- // Add the 2 triangles using vertices 0,1,2 and 2,1,3 (inward normal)
- fprintf(funv, "%10d%s\n", ++nbTria + nbTetra, sELT_TRIA3_DESC ) ;
- fprintf(funv, "%10d%10d%10d\n", i_vert[0]+1, i_vert[1]+1, i_vert[2]+1);
- fprintf(funv, "%10d%s\n", ++nbTria + nbTetra, sELT_TRIA3_DESC ) ;
- fprintf(funv, "%10d%10d%10d\n", i_vert[2]+1, i_vert[1]+1, i_vert[3]+1);
+
+ // Add the 2 triangles using vertices defined in defTria
+ for (int i_tria=0 ; i_tria<2 ; i_tria++) {
+ int i_vert_t[3]; // vertices of the current triangle, in [0..n_z*n_l[
+ for (int i=0 ; i<3 ; i++)
+ i_vert_t[i] = i_vert[defTria[typeHexa][i_tria][i]];
+ fprintf(funv, "%10d%s\n", ++nbTria + nbTetra, sELT_TRIA3_DESC ) ;
+ fprintf(funv, "%10d%10d%10d\n", i_vert_t[0]+1, i_vert_t[1]+1, i_vert_t[2]+1);
+ }
}
fprintf(funv, "%s\n",sUNV_SEPARATOR ) ;
/* ------ */
@@ -334,7 +355,7 @@ int createFemMatrix() {
Mpf_printf(MPI_COMM_WORLD, "Done.\n") ;
}
- Mpf_printf(MPI_COMM_WORLD, "Number of hexaedres = %d\n", nbHexa);
+ Mpf_printf(MPI_COMM_WORLD, "Number of Tetra = %d\n", nbTetra);
Mpf_printf(MPI_COMM_WORLD, "Number of NNZ before sort = %zu\n", nnz_FEM);
// Passer en multithread ?
ierr=Z_sortAndCompactL(z_matComp_FEM, &nnz_FEM) ; CHKERRQ(ierr) ;
diff --git a/src/prepareTEST.c b/src/prepareTEST.c
index 4c7950d..4344511 100644
--- a/src/prepareTEST.c
+++ b/src/prepareTEST.c
@@ -13,6 +13,7 @@ int prepareTEST(void) {
Mpf_printf(MPI_COMM_WORLD, "<PERFTESTS> StepMesh = %e\n" , step);
Mpf_printf(MPI_COMM_WORLD, "<PERFTESTS> NbPts = %d \n", nbPts);
Mpf_printf(MPI_COMM_WORLD, "<PERFTESTS> NbPtsBEM = %d \n", nbPtsBEM);
+ Mpf_printf(MPI_COMM_WORLD, "<PERFTESTS> NbPtsFEM = %d \n", nbPts-nbPtsBEM);
Mpf_printf(MPI_COMM_WORLD, "<PERFTESTS> NbRhs = %d \n", nbRHS);
Mpf_printf(MPI_COMM_WORLD, "<PERFTESTS> nbPtLambda = %e \n", ptsPerLambda);
if (useDetail)
@@ -48,4 +49,3 @@ int displayArray(char *s, void *f, int m, int n) {
Mpf_printf(MPI_COMM_WORLD, "%s[%d, %d]=%e\n", s, i, j, d_f[(size_t)j*m+i]) ;
return 0 ;
}
-
diff --git a/src/rhs.c b/src/rhs.c
index 7b9716a..59831a7 100644
--- a/src/rhs.c
+++ b/src/rhs.c
@@ -1,8 +1,9 @@
#include "main.h"
+#include <mpf_simd.h>
extern double height;
-int computeRhs(void) {
+MPF_SIMD_DISP(int computeRhs(void)) {
double coord[3], OP[3] ;
int i, j, ierr ;
double complex *z_rhs=rhs ;
@@ -21,7 +22,7 @@ int computeRhs(void) {
// OP = fake incident plane wave coming from direction rotating around the object
ierr = computeCoord( (size_t)j*nbPts/nbRHS%nbPts , &(OP[0])) ; CHKERRQ(ierr);
OP[2] -= height/2.; // to have plane waves coming from below too
- z_rhs[(size_t)j * nbPts + i] = cexp(
+ z_rhs[(size_t)j * nbPts + i] = mpf_simd_cexp(
I * k * (coord[0] * OP[0] + coord[1] * OP[1] + coord[2] * OP[2]));
} else {
d_rhs[(size_t)j*nbPts+i]=cos(coord[0]/(j+1)+coord[1]+sin(j)*coord[2])+sin(j) ;
diff --git a/src/testHMAT.c b/src/testHMAT.c
index b3ff1d9..86e9fda 100644
--- a/src/testHMAT.c
+++ b/src/testHMAT.c
@@ -386,10 +386,10 @@ _refineLoop: // Iterative Refinement loop
break;
} /* switch(testedAlgo) */
- /* Destroys the HMAT structure */
#ifdef TESTFEMBEM_DUMP_HMAT
hi->dump_info(hmatrix, "testHMAT_matrix");
#endif
+ /* Destroys the HMAT structure */
HMAT_destroy_matrix( hi, hdesc );
MpfFree(solCLA);
hmat_tracing_dump("testHMAT_trace.json");
--
GitLab