HXT 3D mesher included in gmsh

parent f6300611
......@@ -1020,6 +1020,8 @@ void GetOptions(int argc, char *argv[], bool readConfigFiles, bool exitOnError)
CTX::instance()->mesh.algo2d = ALGO_2D_BAMG;
else if(!strncmp(argv[i], "del3d", 5) || !strncmp(argv[i], "gmsh3d", 6))
CTX::instance()->mesh.algo3d = ALGO_3D_DELAUNAY;
else if(!strncmp(argv[i], "hxt", 3) )
CTX::instance()->mesh.algo3d = ALGO_3D_HXT;
else if(!strncmp(argv[i], "front3d", 7) || !strncmp(argv[i], "netgen", 6))
CTX::instance()->mesh.algo3d = ALGO_3D_FRONTAL;
else if(!strncmp(argv[i], "mmg3d", 5))
......
......@@ -247,6 +247,7 @@
#define ALGO_3D_FRONTAL_HEX 6
#define ALGO_3D_MMG3D 7
#define ALGO_3D_RTREE 9
#define ALGO_3D_HXT 10
// Meshing methods
#define MESH_NONE 0
......
......@@ -15,6 +15,7 @@ set(SRC
meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp
meshGFaceLloyd.cpp meshGFaceOptimize.cpp
meshGRegion.cpp
meshGRegionHxt.cpp
meshDiscreteRegion.cpp
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
......
......@@ -1964,7 +1964,7 @@ static bool buildConsecutiveListOfVertices(
return true;
}
static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
static bool meshGeneratorPeriodic(GFace *gf, bool repairSelfIntersecting1dMesh, bool debug = true)
{
if(CTX::instance()->debugSurface > 0 &&
gf->tag() != CTX::instance()->debugSurface) {
......@@ -2382,6 +2382,9 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
e->g = &CLASS_E;
}
// std::vector<EdgeToRecover> edgesNotRecovered;
for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++) {
std::vector<BDS_Point *> &edgeLoop_BDS = edgeLoops_BDS[i];
for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++) {
......@@ -2389,6 +2392,9 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
edgeLoop_BDS[j]->iD, edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD,
_fatallyFailed);
if(!e) {
// edgesNotRecovered.push_back(EdgeToRecover(edgeLoop_BDS[j]->iD,
// edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD));
Msg::Error("Impossible to recover the edge %d %d", edgeLoop_BDS[j]->iD,
edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
gf->meshStatistics.status = GFace::FAILED;
......@@ -2847,7 +2853,7 @@ void meshGFace::operator()(GFace *gf, bool print)
if(gf->getNativeType() != GEntity::GmshModel &&
(gf->periodic(0) || gf->periodic(1) || singularEdges)) {
if(!meshGeneratorPeriodic(gf, debugSurface >= 0 || debugSurface == -100))
if(!meshGeneratorPeriodic(gf, repairSelfIntersecting1dMesh, debugSurface >= 0 || debugSurface == -100))
Msg::Error("Impossible to mesh periodic face %d", gf->tag());
}
else {
......
......@@ -8,6 +8,7 @@
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "meshGRegion.h"
#include "meshGRegionHxt.h"
#include "meshGFace.h"
#include "meshGFaceOptimize.h"
#include "boundaryLayersData.h"
......@@ -306,7 +307,12 @@ void MeshDelaunayVolume(std::vector<GRegion *> &regions)
// now do insertion of points
if(CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D) {
if(CTX::instance()->mesh.algo3d == ALGO_3D_HXT) {
if (meshGRegionHxt (gr) != 0){
Msg::Error ("HXT 3D mesh failed");
}
}
else if(CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D) {
refineMeshMMG(gr);
}
else if(CTX::instance()->mesh.oldRefinement) {
......
......@@ -8,15 +8,32 @@ set(SRC
hxt_context.c
hxt_curvature.c
hxt_edge.c
hxt_laplace.c
hxt_linear_system.c
hxt_linear_system_lu.c
hxt_linear_system_petsc.c
hxt_mesh.c
hxt_message.c
hxt_non_linear_solver.c
hxt_option.c
hxt_parametrization.c
hxt_mean_values.c
hxt_opt.c
hxt_sort.c
hxt_tools.c
hxt_mesh.c
predicates.c
hxt_mesh3d.c
hxt_mesh3d_main.c
hxt_mesh_size.c
hxt_tetOpti.c
hxt_tetrahedra.c
hxt_tetRepair.c
hxt_tetUtils.c
hxt_tetFlag.c
hxt_tetPostpro.c
hxt_tet_aspect_ratio.c
hxt_vertices.c
hxt_parametrization.c
hxt_mean_values.c
hxt_boundary_recovery.cxx
)
file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
......
Hxt copyright (C) 2017-2018
Universite catholique de Louvain
The TetGen/BR code (tetgenBR.{cxx,h}) is copyright (c) 2016 Hang Si,
Weierstrass Institute for Applied Analysis and Stochatics. It is relicensed
under the terms of gmsh/LICENSE.txt for use in Gmsh thanks to a Software License
Agreement between Weierstrass Institute for Applied Analysis and Stochastics and
GMESH SPRL.
This diff is collapsed.
......@@ -29,12 +29,17 @@ typedef enum
HXT_STATUS_OUT_OF_MEMORY = -4,
HXT_STATUS_FILE_CANNOT_BE_OPENED = -5,
HXT_STATUS_POINTER_ERROR = -6,
HXT_STATUS_READ_ERROR = -7,
HXT_STATUS_WRITE_ERROR = -8,
HXT_STATUS_RANGE_ERROR = -9,
HXT_STATUS_FORMAT_ERROR = -10,
// INTERNAL Errors (<= HXT_STATUS_INTERNAL) => HXT_CHECK does not give trace message but returns... should be catched internally !
HXT_STATUS_INTERNAL = -1024,
HXT_STATUS_SKIP = -1025,
HXT_STATUS_TRYAGAIN = -1026
HXT_STATUS_CONFLICT = -1025,
HXT_STATUS_SKIP = -1026,
HXT_STATUS_TRYAGAIN = -1027
}HXTStatus;
......
#include "hxt_bbox.h"
// /* create a new bounding box */
// HXTStatus hxtBboxCreate(hxtBbox** bboxP){
// HXTStatus hxtBboxCreate(HXTBbox** bboxP){
// HXT_CHECK(
// hxtMalloc((void **) bboxP, sizeof(hxtBbox)) );
// hxtMalloc((void **) bboxP, sizeof(HXTBbox)) );
// unsigned i;
// for (i=0; i<3; i++)
......@@ -16,7 +16,7 @@
// }
// HXTStatus hxtBboxDelete(hxtBbox** bboxP){
// HXTStatus hxtBboxDelete(HXTBbox** bboxP){
// HXT_CHECK(
// hxtFree((void **) bboxP) );
......@@ -25,7 +25,7 @@
/* update the bounding box with one new vertex */
HXTStatus hxtBboxAddOne(hxtBbox* bbox, double* coord){
HXTStatus hxtBboxAddOne(HXTBbox* bbox, double* coord){
unsigned i;
for (i=0; i<3; i++)
{
......@@ -41,7 +41,7 @@ HXTStatus hxtBboxAddOne(hxtBbox* bbox, double* coord){
/* update the bounding box with an array of n vertices at once (far quicker) */
HXTStatus hxtBboxAdd(hxtBbox* bbox, double* coord, const uint32_t n){
HXTStatus hxtBboxAdd(HXTBbox* bbox, double* coord, const uint32_t n){
#pragma omp parallel
{
......@@ -51,19 +51,19 @@ HXTStatus hxtBboxAdd(hxtBbox* bbox, double* coord, const uint32_t n){
#pragma omp for nowait
for (uint32_t i=0; i<n; i++)
{
if(coord[4*i ]<min[0])
min[0]=coord[4*i ];
if(coord[4*i+1]<min[1])
min[1]=coord[4*i+1];
if(coord[4*i+2]<min[2])
min[2]=coord[4*i+2];
if(coord[4*i ]>max[0])
max[0]=coord[4*i ];
if(coord[4*i+1]>max[1])
max[1]=coord[4*i+1];
if(coord[4*i+2]>max[2])
max[2]=coord[4*i+2];
if(coord[(size_t) 4*i ]<min[0])
min[0]=coord[(size_t) 4*i ];
if(coord[(size_t) 4*i+1]<min[1])
min[1]=coord[(size_t) 4*i+1];
if(coord[(size_t) 4*i+2]<min[2])
min[2]=coord[(size_t) 4*i+2];
if(coord[(size_t) 4*i ]>max[0])
max[0]=coord[(size_t) 4*i ];
if(coord[(size_t) 4*i+1]>max[1])
max[1]=coord[(size_t) 4*i+1];
if(coord[(size_t) 4*i+2]>max[2])
max[2]=coord[(size_t) 4*i+2];
}
#pragma omp critical
......
......@@ -12,15 +12,15 @@ extern "C" {
typedef struct hxtBboxStruct{
double hxtDeclareAligned min[3];
double hxtDeclareAligned max[3];
} hxtBbox;
} HXTBbox;
// /* create a new bounding box */
// HXTStatus hxtBboxCreate(hxtBbox** bboxP);
// HXTStatus hxtBboxCreate(HXTBbox** bboxP);
// /* delete a bounding box */
// HXTStatus hxtBboxDelete(hxtBbox** bboxP);
// HXTStatus hxtBboxDelete(HXTBbox** bboxP);
static inline void hxtBboxInit(hxtBbox* bbox){
static inline void hxtBboxInit(HXTBbox* bbox){
bbox->min[0] = DBL_MAX;
bbox->min[1] = DBL_MAX;
bbox->min[2] = DBL_MAX;
......@@ -30,10 +30,10 @@ static inline void hxtBboxInit(hxtBbox* bbox){
}
/* update the bounding box with one new vertex */
HXTStatus hxtBboxAddOne(hxtBbox* bbox, double* coord);
HXTStatus hxtBboxAddOne(HXTBbox* bbox, double* coord);
/* update the bounding box with an array of n vertices at once (far quicker) */
HXTStatus hxtBboxAdd(hxtBbox* bbox, double* coord, uint32_t n);
HXTStatus hxtBboxAdd(HXTBbox* bbox, double* coord, uint32_t n);
#ifdef __cplusplus
}
......
This diff is collapsed.
#ifndef _HXT_BOUNDARY_RECOVERY_
#define _HXT_BOUNDARY_RECOVERY_
HXTStatus hxt_boundary_recovery(HXTMesh *mesh);
#endif
......@@ -2,6 +2,7 @@
#define HEXTREME_LAPLACE_H
#include "hxt_mesh.h"
HXTStatus hxtLaplace(HXTMesh *mesh);
#endif
......@@ -114,24 +114,34 @@ typedef hxtDeclareAligned double alignedDouble;
// y[from:to] += a*x[from:to]
// y and x must be 256-bit aligned
// memory should be allocated in consequence
static void rowAxpy(double a, alignedDouble *__restrict__ x, alignedDouble *__restrict__ y, int from, int to)
static void rowAxpy(double a, double *__restrict__ px, double *__restrict__ py, int from, int to)
{
// Can't use the aligned attribute on function arguments.
hxtDeclareAligned double * __restrict__ x = px;
hxtDeclareAligned double * __restrict__ y = py;
int i = from;
int pfrom = (from+3)&(~3);
int pfrom = (from+7)&(~7);
if (pfrom > to)
pfrom = to;
for (; i < pfrom; ++i)
y[i] += a*x[i];
for (; i+15 < to; i+=16) {
alignedDouble * __restrict__ X = x + i;
alignedDouble * __restrict__ Y = y + i;
hxtDeclareAligned double * __restrict__ X = x + i;
hxtDeclareAligned double * __restrict__ Y = y + i;
for (int j = 0; j < 16; ++j){
Y[j] += a * X[j];
}
}
for (; i+7 < to; i+=8) {
hxtDeclareAligned double * __restrict__ X = x + i;
hxtDeclareAligned double * __restrict__ Y = y + i;
for (int j = 0; j < 8; ++j)
Y[j] += a * X[j];
}
for (; i+3 < to; i+=4) {
alignedDouble * __restrict__ X = x + i;
alignedDouble * __restrict__ Y = y + i;
double * __restrict__ X = x + i;
double * __restrict__ Y = y + i;
for (int j = 0; j < 4; ++j)
Y[j] += a * X[j];
}
......@@ -144,35 +154,40 @@ static int imin(int a, int b) {
return a < b ? a : b;
}
static double rowReduction(alignedDouble *__restrict__ x, alignedDouble *__restrict__ y, int from, int to)
static double rowReduction(double *__restrict__ px, double *__restrict__ py, int from, int to)
{
int i = from;
double r = 0;
int pfrom = (from+3)&(~3);
hxtDeclareAligned double * __restrict__ x = px;
hxtDeclareAligned double * __restrict__ y = py;
int i = from;
int pfrom = (from+7)&(~7);
for (; i < imin(pfrom, to); ++i)
r += x[i] * y[i];
double R[4];
for (; i+3 < to; i+=4) {
alignedDouble * __restrict__ X = x + i;
alignedDouble * __restrict__ Y = y + i;
for (int j = 0; j < 4; ++j)
double R[8];
for (; i+7 < to; i+=8) {
hxtDeclareAligned double * __restrict__ X = x + i;
hxtDeclareAligned double * __restrict__ Y = y + i;
for (int j = 0; j < 8; ++j)
R[j] = X[j]*Y[j];
r += R[0]+R[1]+R[2]+R[3];
r += R[0]+R[1]+R[2]+R[3]+R[4]+R[5]+R[6]+R[7];
}
for (; i < to; ++i)
r += x[i] * y[i];
return r;
}
static void rowZero(alignedDouble *__restrict__ x, int from, int to)
static void rowZero(double *__restrict__ px, int from, int to)
{
hxtDeclareAligned double * __restrict__ x = px;
int i = from;
int pfrom = (from+3)&(~3);
int pfrom = (from+7)&(~7);
for (; i < imin(pfrom, to); ++i)
x[i] = 0;
for (; i+3 < to; i+=4) {
alignedDouble * __restrict__ X = x + i;
for (int j = 0; j < 4; ++j)
for (; i+7 < to; i+=8) {
hxtDeclareAligned double * __restrict__ X = x + i;
for (int j = 0; j < 8; ++j)
X[j] = 0;
}
for (; i < to; ++i)
......@@ -279,7 +294,7 @@ HXTStatus hxtLinearSystemLUCreate(HXTLinearSystemLU **pSystem, int nElements, in
}
free(nodeRowStart);
free(nodeRowEnd);
system->M = malloc(sizeof(double)*totalSize); // FIXME Gmsh instead of _mm_malloc
system->M = _mm_malloc(sizeof(double)*totalSize, PADDING*8);
system->rows = malloc(sizeof(double*)*system->n);
for (int i = 0; i < totalSize; ++i)
system->M[i] = 0;
......@@ -291,7 +306,7 @@ HXTStatus hxtLinearSystemLUCreate(HXTLinearSystemLU **pSystem, int nElements, in
totalSize += system->rowEnd[i]-system->rowStart[i]+(paddedStart-start);
system->rows[i] = system->M + paddedStart;
}
system->x = malloc(sizeof(double)*system->n); // FIXME Gmsh instead of _mm_malloc
system->x = _mm_malloc(sizeof(double)*system->n, PADDING*8);
return HXT_STATUS_OK;
}
......@@ -384,7 +399,7 @@ HXTStatus hxtLinearSystemLUAddMatrixEntry(HXTLinearSystemLU *system, int node0,
HXT_ERROR_MSG(HXT_STATUS_FAILED, "node %i or %i not in the domain", node0, node1);
int row0 = system->nodeMap[node0]*system->nFields + field0;
int col1 = system->nodeMap[node1]*system->nFields + field1;
system->rows[row0][col1] += v;
return HXT_STATUS_OK;
}
......@@ -403,7 +418,7 @@ HXTStatus hxtLinearSystemLUSolve(HXTLinearSystemLU *system, double *rhs, double
LUPDecompose(system);
system->flaglu=1;
}
LUPSolve(system, rhs);
for (int i = 0; i < system->nNodes; ++i){
int ii = system->nodeMap[i];
......
......@@ -229,16 +229,30 @@ HXTStatus hxtLinearSystemPETScAddRhsEntry(HXTLinearSystemPETSc *system, double *
HXTStatus hxtLinearSystemPETScSolve(HXTLinearSystemPETSc *system, double *rhs, double *solution){
Vec b;
Vec x;
HXT_PETSC_CHECK(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, system->nDofs, rhs, &b));
HXT_PETSC_CHECK(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, system->nDofs, solution, &x));
double normTest=0.0;
for (int iv = 0; iv < system->nDofs; ++iv){
normTest+=solution[iv]*solution[iv];
}
/* printf("IN PETSC norm vect sol : %g ; nofs : %i\n",normTest,system->nDofs); */
if(system->assemblyNeeded) {
HXT_PETSC_CHECK(MatAssemblyBegin(system->a, MAT_FINAL_ASSEMBLY));
HXT_PETSC_CHECK(MatAssemblyEnd(system->a, MAT_FINAL_ASSEMBLY));
system->assemblyNeeded = 0;
}
HXT_PETSC_CHECK(KSPSetOperators(system->ksp, system->a, system->a));
HXT_PETSC_CHECK(KSPSolve(system->ksp, b, system->x));
/* HXT_PETSC_CHECK(KSPSolve(system->ksp, b, system->x)); */
KSPSetInitialGuessNonzero(system->ksp,PETSC_TRUE);
PetscBool flag[1];
KSPGetInitialGuessNonzero(system->ksp,flag);
/* printf("IN PETSC initialguessnonzero : %i\n",flag[0]); */
HXT_PETSC_CHECK(KSPSolve(system->ksp, b, x));
HXT_PETSC_CHECK(VecDestroy(&b));
HXT_CHECK(hxtLinearSystemPETScMapFromVec(system, system->x, solution));
HXT_CHECK(hxtLinearSystemPETScMapFromVec(system, x, solution));
HXT_PETSC_CHECK(VecDestroy(&x));
/* HXT_CHECK(hxtLinearSystemPETScMapFromVec(system, system->x, solution)); */
return HXT_STATUS_OK;
}
......
......@@ -155,12 +155,12 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
// gather local edge index
for(int it=0; it<2; it++){
if(tri[it]==(uint64_t)-1)
break;
break;
for(int k=0; k<3; k++){
if(edges->tri2edg[3*tri[it]+k]==ie){
ik[it] = k;
break;
}
if(edges->tri2edg[3*tri[it]+k]==ie){
ik[it] = k;
break;
}
}
}
......@@ -170,32 +170,32 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
uint32_t v1 = edges->node[2*ie+(1-ij)];
if (flag[v0]==1){//boundary nodes/conditons
HXT_CHECK(hxtLinearSystemSetMatrixRowIdentity(sys,v0,0));
HXT_CHECK(hxtLinearSystemSetRhsEntry(sys,Urhs, v0,0, uv[2*v0+0]));
HXT_CHECK(hxtLinearSystemSetRhsEntry(sys,Vrhs, v0,0, uv[2*v0+1]));
HXT_CHECK(hxtLinearSystemSetMatrixRowIdentity(sys,v0,0));
HXT_CHECK(hxtLinearSystemSetRhsEntry(sys,Urhs, v0,0, uv[2*v0+0]));
HXT_CHECK(hxtLinearSystemSetRhsEntry(sys,Vrhs, v0,0, uv[2*v0+1]));
}
else {// inner node
double e[3] = {mesh->vertices.coord[4*v1+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*v1+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*v1+2]-mesh->vertices.coord[4*v0+2]};
double ne = sqrt(e[0]*e[0]+e[1]*e[1]+e[2]*e[2]);
uint32_t vLeft = mesh->triangles.node[3*tri[0] + (ik[0]+2)%3];
double a[3] = {mesh->vertices.coord[4*vLeft+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*vLeft+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*vLeft+2]-mesh->vertices.coord[4*v0+2]};
double na = sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
double thetaL = acos((a[0]*e[0]+a[1]*e[1]+a[2]*e[2])/(na*ne));
double thetaR=0.;
if(tri[1]==(uint64_t)-1)
thetaR=0.;
else{
uint32_t vRight = mesh->triangles.node[3*tri[1] + (ik[1]+2)%3];
double b[3] = {mesh->vertices.coord[4*vRight+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*vRight+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*vRight+2]-mesh->vertices.coord[4*v0+2]};
double nb = sqrt(b[0]*b[0]+b[1]*b[1]+b[2]*b[2]);
thetaR = acos((b[0]*e[0]+b[1]*e[1]+b[2]*e[2])/(nb*ne));
}
double c = (tan(.5*thetaL)+tan(.5*thetaR))/ne;
HXT_CHECK(hxtLinearSystemAddMatrixEntry(sys, v0, 0, v1, 0, -c));
HXT_CHECK(hxtLinearSystemAddMatrixEntry(sys, v0, 0, v0, 0, c));
else {// inner node
double e[3] = {mesh->vertices.coord[4*v1+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*v1+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*v1+2]-mesh->vertices.coord[4*v0+2]};
double ne = sqrt(e[0]*e[0]+e[1]*e[1]+e[2]*e[2]);
uint32_t vLeft = mesh->triangles.node[3*tri[0] + (ik[0]+2)%3];
double a[3] = {mesh->vertices.coord[4*vLeft+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*vLeft+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*vLeft+2]-mesh->vertices.coord[4*v0+2]};
double na = sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
double thetaL = acos((a[0]*e[0]+a[1]*e[1]+a[2]*e[2])/(na*ne));
double thetaR=0.;
if(tri[1]==(uint64_t)-1)
thetaR=0.;
else{
uint32_t vRight = mesh->triangles.node[3*tri[1] + (ik[1]+2)%3];
double b[3] = {mesh->vertices.coord[4*vRight+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*vRight+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*vRight+2]-mesh->vertices.coord[4*v0+2]};
double nb = sqrt(b[0]*b[0]+b[1]*b[1]+b[2]*b[2]);
thetaR = acos((b[0]*e[0]+b[1]*e[1]+b[2]*e[2])/(nb*ne));
}
double c = (tan(.5*thetaL)+tan(.5*thetaR))/ne;
HXT_CHECK(hxtLinearSystemAddMatrixEntry(sys, v0, 0, v1, 0, -c));
HXT_CHECK(hxtLinearSystemAddMatrixEntry(sys, v0, 0, v0, 0, c));
}
}
}// end for int ie
......@@ -224,7 +224,7 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
HXT_CHECK(hxtLinearSystemAddMatrixEntry(sys, edges->node[2*edges_ll[j]+0], 0, edges->node[2*edges_ll[j]+1], 0, -tan(.5*(M_PI-theta)/2)/le ));
HXT_CHECK(hxtLinearSystemAddMatrixEntry(sys, edges->node[2*edges_ll[j]+1], 0, edges->node[2*edges_ll[j]+1], 0, tan(.5*(M_PI-theta)/2)/le ));
HXT_CHECK(hxtLinearSystemAddMatrixEntry(sys, edges->node[2*edges_ll[j]+1], 0, edges->node[2*edges_ll[j]+0], 0, -tan(.5*(M_PI-theta)/2)/le ));
c[j] = tan(.5*(M_PI-theta)/2)/radius;
d[j] = tan(.5*(theta))/radius;
}
......@@ -243,7 +243,7 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
}
for(int j=0; j<n; j++){
uint32_t v = edges->node[2*edges_ll[j]+0];
HXT_CHECK(hxtLinearSystemAddMatrixEntry(sys, v, 0, v, 0, c[j]));
......@@ -281,22 +281,60 @@ HXTStatus hxtMeanValueAspectRatio(HXTMeanValues *meanValues, int *aspectRatio)
{
if(meanValues->aspectRatio<0){
*aspectRatio = 1;
HXTMesh *mesh = meanValues->initialEdges->edg2mesh;
double grad[3][2] = {{-1./2.,-sqrt(3)/6.},{1./2.,-sqrt(3)/6.},{0.,sqrt(3)/3.}};
uint64_t nTri = mesh->triangles.num;
double *uv = meanValues->uv;
for(uint64_t it=0; it<nTri; it++){
uint32_t *nodes = mesh->triangles.node + 3*it;
double jac[2][2] = {{0.,0.},{0.,0.}};
for(int i=0; i<3; i++){
jac[0][0] += uv[2*nodes[i]+0]*grad[i][0];// dx dxi
jac[0][1] += uv[2*nodes[i]+0]*grad[i][1];// dx deta
jac[1][0] += uv[2*nodes[i]+1]*grad[i][0];// dy dxi
jac[1][1] += uv[2*nodes[i]+1]*grad[i][1];// dy deta
}
double det = jac[0][0]*jac[1][1]-jac[1][0]*jac[0][1];
double frob = 0.;
for(int i =0; i<2; i++)
for(int j=0; j<2; j++)
frob += jac[i][j]*jac[i][j];
double quality = 2*det/frob;
if(quality<.1){
printf("wrong aspect ratio !!!!!!! D-: \t %f\n",quality);
*aspectRatio = 0;
break;
}
}
/*
*aspectRatio = 1;
uint32_t numEdges = meanValues->initialEdges->numEdges;
double *uv = meanValues->uv;
uint32_t *nodes = meanValues->initialEdges->node;
for(uint32_t i=0; i<numEdges; i++){
double du = uv[2*nodes[2*i+1]+0] - uv[2*nodes[2*i+0]+0];
double dv = uv[2*nodes[2*i+1]+1] - uv[2*nodes[2*i+0]+1];
if(sqrt(du*du+dv*dv) < 1e-4){
*aspectRatio = 0;
break;
}
}
double du = uv[2*nodes[2*i+1]+0] - uv[2*nodes[2*i+0]+0];
double dv = uv[2*nodes[2*i+1]+1] - uv[2*nodes[2*i+0]+1];
if(sqrt(du*du+dv*dv) < 1e-4){
*aspectRatio = 0;
break;
}
}
*/
}
else
*aspectRatio = meanValues->aspectRatio;
......@@ -329,7 +367,7 @@ HXTStatus hxtMeanValuesGetData(HXTMeanValues *mv, uint64_t **global,uint32_t **g
global_[ie] = mv->initialEdges->global[ie];
if(gn!=NULL)
for(int kk=0; kk<3; kk++)
gn_[3*ie+kk] = mv->initialEdges->edg2mesh->triangles.node[3*ie+kk];
gn_[3*ie+kk] = mv->initialEdges->edg2mesh->triangles.node[3*ie+kk];
}
......
This diff is collapsed.
......@@ -4,13 +4,22 @@
#include "hxt_tools.h" // to have SIMD_ALIGN and stdint.h
#define HXT_GHOST_VERTEX UINT32_MAX
#define HXT_DELETED_COLOR (UINT16_MAX-1)
#define HXT_NO_ADJACENT UINT64_MAX
/* Element types, same as Gmsh */
typedef enum {
HXT_NO_ELT = 0,
HXT_LINE = 1,
HXT_TRI = 2,
HXT_QUAD = 3,
HXT_TET = 4,
HXT_HEX = 5,
HXT_PRI = 6,
HXT_PYR = 7
} HXT_ELT_TYPE;
struct hxtMeshStruct {
HXTContext* ctx;
......@@ -21,25 +30,68 @@ struct hxtMeshStruct {
uint32_t size;
double* coord; // 3 coordinates + 1 double per vertex!
} vertices;
uint64_t* typeStat; // to delete
// tetrahedra
struct {
uint32_t* node; // aligned (size = tetrahedra.size*4*sizeof(uint32_t))
uint64_t* neigh; // aligned (size = tetrahedra.size*4*sizeof(uint64_t))
char* neighType;
uint16_t* colors;
uint16_t* flag;
uint64_t num; // number of tetrahedra
uint64_t size; // reserved number of tetrahedra (size of the vector)
} tetrahedra;
// hexahedra
struct {
uint32_t* node; // aligned (size = hexahedra.size*8*sizeof(uint32_t))
uint64_t* neigh; // aligned (size = hexahedra.size*6*sizeof(uint64_t))
char* neighType;
uint