hxt 3D mesher properly connected

parent 18ff77e9
This diff is collapsed.
......@@ -24,17 +24,31 @@ set(SRC
hxt_mesh3d_main.c
hxt_mesh_size.c
hxt_tetOpti.c
hxt_tetrahedra.c
hxt_tetRepair.c
hxt_tetUtils.c
hxt_tetColor.c
hxt_tetFlag.c
hxt_tetPostpro.c
hxt_tet_aspect_ratio.c
hxt_tetDelaunay.c
hxt_vertices.c
hxt_parametrization.c
hxt_mean_values.c
hxt_boundary_recovery.cxx
)
# do not use arithmetic contraction in predicates.c
if ("x${CMAKE_C_COMPILER_ID}" STREQUAL "xMSVC" )
set_source_files_properties(predicates.c PROPERTIES COMPILE_FLAGS "/fp:strict")
endif()
if (CMAKE_C_COMPILER_ID MATCHES "GNU|Clang")
set_source_files_properties(predicates.c PROPERTIES COMPILE_FLAGS "-fno-unsafe-math-optimizations -ffp-contract=off")
endif()
if (CMAKE_C_COMPILER_ID STREQUAL "Intel")
set_source_files_properties(predicates.c PROPERTIES COMPILE_FLAGS "-fp-model strict")
endif()
file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
append_gmsh_src(contrib/hxt "${SRC};${HDR}")
......@@ -82,3 +82,22 @@ HXTStatus hxtBboxAdd(HXTBbox* bbox, double* coord, const uint32_t n){
return HXT_STATUS_OK;
}
HXTStatus hxtBboxMerge(HXTBbox* bbox1, HXTBbox* bbox2, HXTBbox* bboxResult)
{
unsigned i;
for (i=0; i<3; i++)
{
if(bbox1->min[i]<bbox2->min[i])
bboxResult->min[i] = bbox1->min[i];
else
bboxResult->min[i] = bbox2->min[i];
if(bbox1->max[i]>bbox2->max[i])
bboxResult->max[i] = bbox1->max[i];
else
bboxResult->max[i] = bbox2->max[i];
}
return HXT_STATUS_OK;
}
......@@ -35,6 +35,9 @@ 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);
/* merge two bbox (result can be a pointer to bbox1 or bbox2) */
HXTStatus hxtBboxMerge(HXTBbox* bbox1, HXTBbox* bbox2, HXTBbox* bboxResult);
#ifdef __cplusplus
}
#endif
......
This diff is collapsed.
......@@ -294,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;
......@@ -306,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;
}
......@@ -399,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;
}
......@@ -418,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,30 +229,16 @@ 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)); */
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(KSPSolve(system->ksp, b, system->x));
HXT_PETSC_CHECK(VecDestroy(&b));
HXT_CHECK(hxtLinearSystemPETScMapFromVec(system, x, solution));
HXT_PETSC_CHECK(VecDestroy(&x));
/* HXT_CHECK(hxtLinearSystemPETScMapFromVec(system, system->x, solution)); */
HXT_CHECK(hxtLinearSystemPETScMapFromVec(system, system->x, solution));
return HXT_STATUS_OK;
}
......
......@@ -77,6 +77,18 @@ HXTStatus hxtMeshCreate ( HXTContext* ctx, HXTMesh** mesh) {
(*mesh)->lines.num = 0;
(*mesh)->lines.size = 0;
// boundary representation
(*mesh)->brep.numVolumes = 0;
(*mesh)->brep.numSurfacesPerVolume = NULL;
(*mesh)->brep.surfacesPerVolume = NULL;
(*mesh)->brep.numSurfaces = 0;
(*mesh)->brep.numCurvesPerSurface = NULL;
(*mesh)->brep.curvesPerSurface = NULL;
(*mesh)->brep.numCurves = 0;
(*mesh)->brep.endPointsOfCurves = NULL;
(*mesh)->brep.numPoints = 0;
(*mesh)->brep.points = NULL;
return HXT_STATUS_OK;
}
......@@ -127,6 +139,14 @@ HXTStatus hxtMeshDelete ( HXTMesh** mesh) {
HXT_CHECK( hxtAlignedFree(&(*mesh)->lines.node) );
HXT_CHECK( hxtAlignedFree(&(*mesh)->lines.colors) );
// boundary representation
HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.numSurfacesPerVolume) );
HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.surfacesPerVolume) );
HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.numCurvesPerSurface) );
HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.curvesPerSurface) );
HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.endPointsOfCurves) );
HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.points) );
HXT_CHECK( hxtFree(mesh) );
return HXT_STATUS_OK;
......@@ -450,7 +470,7 @@ HXTStatus hxtMeshWriteGmsh ( HXTMesh* mesh , const char *filename) {
if(mesh->tetrahedra.node[i*4 + 3]!=UINT32_MAX){
uint16_t myColor = mesh->tetrahedra.colors ? mesh->tetrahedra.colors[i] : 0;
// color = UINT16_MAX --> outside the domain
if (myColor != UINT16_MAX)
// if (myColor != UINT16_MAX)
++index;
}
}
......@@ -498,7 +518,7 @@ HXTStatus hxtMeshWriteGmsh ( HXTMesh* mesh , const char *filename) {
for (i=0; i<mesh->tetrahedra.num; i++){
if(mesh->tetrahedra.node[i*4 + 3]!=UINT32_MAX){
uint16_t myColor = mesh->tetrahedra.colors ? mesh->tetrahedra.colors[i] : 0;
if (myColor != UINT16_MAX)
// if (myColor != UINT16_MAX)
fprintf(file,"%lu %u 2 0 %u %u %u %u %u\n", ++index,TETID,
myColor,
mesh->tetrahedra.node[i*4]+1,
......
......@@ -3,8 +3,8 @@
#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
......@@ -34,18 +34,18 @@ struct hxtMeshStruct {
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;
uint8_t* 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;
uint8_t* neighType;
uint16_t* colors;
uint16_t* flag;
uint64_t num; // number of tetrahedra
......@@ -56,18 +56,18 @@ struct hxtMeshStruct {
struct {
uint32_t* node; // aligned (size = prisms.size*6*sizeof(uint32_t))
uint64_t* neigh; // aligned (size = prisms.size*5*sizeof(uint64_t))
char* neighType;
uint8_t* neighType;
uint16_t* colors;
uint16_t* flag;
uint64_t num; // number of tetrahedra
uint64_t size; // reserved number of prisms (size of the vector)
} prisms;
// pyramids
struct {
uint32_t* node; // aligned (size = pyramids.size*5*sizeof(uint32_t))
uint64_t* neigh; // aligned (size = pyramids.size*5*sizeof(uint64_t))
char* neighType;
uint8_t* neighType;
uint16_t* colors;
uint16_t* flag;
uint64_t num; // number of tetrahedra
......@@ -77,7 +77,7 @@ struct hxtMeshStruct {
// triangles // TODO: consider writing a array of structure...
struct {
uint32_t* node;
uint64_t* neigh;
uint64_t* neigh;
uint16_t* colors;
uint64_t num;
uint64_t size;
......@@ -98,17 +98,19 @@ struct hxtMeshStruct {
uint64_t num;
uint64_t size;
} lines;
// boundary representation
struct {
uint32_t numVolumes;
uint32_t *numSurfacesPerVolume;
uint32_t *surfacesPerVolume;
uint32_t numSurfaces;
uint32_t *numCurvesPerSurface;
uint32_t *curvesPerSurface;
uint32_t numCurves;
uint32_t *endPointsOfCurves;
uint32_t numPoints;
uint32_t *points;
uint16_t numVolumes;
uint16_t *numSurfacesPerVolume;
uint16_t *surfacesPerVolume;
uint16_t numSurfaces;
uint16_t *numCurvesPerSurface;
uint16_t *curvesPerSurface;
uint16_t numCurves;
uint16_t *endPointsOfCurves;
uint16_t numPoints;
uint16_t *points;
} brep;
};
......
// #include "hxt_mesh_size.h"
#include "hxt_tetrahedra.h"
#include "hxt_tetDelaunay.h"
// #include "hxt_vertices.h"
#include "hxt_mesh3d.h"
#include "predicates.h"
......@@ -17,7 +17,8 @@
HXTStatus hxtComputeMeshSizeFromTrianglesAndLines(HXTMesh* mesh, HXTDelaunayOptions* delOptions) {
HXTStatus hxtCreateNodalsizeFromTrianglesAndLines(HXTMesh* mesh, HXTDelaunayOptions* delOptions)
{
HXT_CHECK(hxtAlignedMalloc(&delOptions->nodalSizes,mesh->vertices.num*sizeof(double)));
......@@ -63,7 +64,8 @@ HXTStatus hxtComputeMeshSizeFromTrianglesAndLines(HXTMesh* mesh, HXTDelaunayOpti
return HXT_STATUS_OK;
}
HXTStatus hxtComputeMeshSizeFromMesh (HXTMesh* mesh, HXTDelaunayOptions* delOptions) {
HXTStatus hxtCreateNodalsizeFromMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions)
{
HXT_CHECK(hxtAlignedMalloc(&delOptions->nodalSizes,mesh->vertices.num*sizeof(double)));
......@@ -92,13 +94,20 @@ HXTStatus hxtComputeMeshSizeFromMesh (HXTMesh* mesh, HXTDelaunayOptions* delOpti
}
}
}
return HXT_STATUS_OK;
return HXT_STATUS_OK;
}
HXTStatus hxtDestroyNodalsize(HXTDelaunayOptions* delOptions)
{
HXT_CHECK( hxtAlignedFree(&delOptions->nodalSizes) );
return HXT_STATUS_OK;
}
HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions){
HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions)
{
// we assume that the input is a surface mesh
if (mesh->tetrahedra.num)
return HXT_ERROR_MSG(HXT_STATUS_FAILED, "The input mesh should only contain triangles");
......@@ -156,78 +165,12 @@ HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions){
/***************************************************
* Coloring the mesh *
***************************************************/
HXTStatus hxtColorMesh(HXTMesh* mesh, uint16_t *nbColors) {
uint64_t *stack;
HXT_CHECK(hxtMalloc(&stack,mesh->tetrahedra.num*sizeof(uint64_t)));
// now that tetrahedra are flaged, we can proceed to colorize the mesh
memset(mesh->tetrahedra.colors, 0, mesh->tetrahedra.size*sizeof(uint16_t));
uint16_t color = 1;
uint16_t colorOut = 0;
while (1){
uint64_t stackSize = 0;
uint64_t first = UINT64_MAX;
for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
if(mesh->tetrahedra.colors[i]==0){
first = i;
break;
}
}
if(first==UINT64_MAX)
break;
stack[stackSize++] = first;
mesh->tetrahedra.colors[first] = color;
for (uint64_t i=0; i<stackSize; i++) {
uint64_t t = stack[i];
if (mesh->tetrahedra.node[4*t+3] == HXT_GHOST_VERTEX)
colorOut = color;
for (unsigned j=0; j<4; j++) {
if(mesh->tetrahedra.neigh[4*t+j]!=HXT_NO_ADJACENT && isFacetConstrained(mesh, 4*t+j)==0){ // the facet is not a boundary facet
uint64_t neigh = mesh->tetrahedra.neigh[4*t+j]/4;
if(mesh->tetrahedra.colors[neigh]==0){
stack[stackSize++] = neigh;
mesh->tetrahedra.colors[neigh] = color;
}
}
}
}
color++;
}
*nbColors = color-1; // -1 because we began at one AND the colorout is counted...
HXT_CHECK( hxtFree(&stack) );
#pragma omp parallel for
for (int i=0;i<mesh->tetrahedra.num;i++){
if (mesh->tetrahedra.colors[i] == colorOut){
mesh->tetrahedra.colors[i] = UINT16_MAX;
}
else if(mesh->tetrahedra.colors[i] > colorOut){
mesh->tetrahedra.colors[i]--;
}
}
return HXT_STATUS_OK;
}
// refine
double hxtTetCircumcenter(double a[3], double b[3], double c[3], double d[3],
double circumcenter[3], double *xi, double *eta, double *zeta){
double circumcenter[3], double *xi, double *eta, double *zeta)
{
double xba, yba, zba, xca, yca, zca, xda, yda, zda;
double balength, calength, dalength;
double xcrosscd, ycrosscd, zcrosscd;
......@@ -326,7 +269,7 @@ HXTStatus hxtRefineTetrahedraOneStep(HXTMesh* mesh, HXTDelaunayOptions* delOptio
for (uint64_t i=0; i<mesh->tetrahedra.num; i++)
{
newVertices[(size_t) 4*i+3] = -1.0;
if (mesh->tetrahedra.colors[i] != UINT16_MAX && isTetProcessed(mesh, i)==0){
if (mesh->tetrahedra.colors[i] != UINT16_MAX && getProcessedFlag(mesh, i)==0){
double *a = mesh->vertices.coord + (size_t) 4*mesh->tetrahedra.node[4*i+0];
double *b = mesh->vertices.coord + (size_t) 4*mesh->tetrahedra.node[4*i+1];
double *c = mesh->vertices.coord + (size_t) 4*mesh->tetrahedra.node[4*i+2];
......@@ -372,7 +315,7 @@ HXTStatus hxtRefineTetrahedraOneStep(HXTMesh* mesh, HXTDelaunayOptions* delOptio
localAdd++;
}
markTetAsProcessed(mesh, i); // we do not need to refine that tetrahedra anymore
setProcessedFlag(mesh, i); // we do not need to refine that tetrahedra anymore
}
}
......
#ifndef _HXT_MESH_3D_
#define _HXT_MESH_3D_
#include "hxt_tetrahedra.h"
#include "hxt_tetDelaunay.h"
/// Creates a structure that allows to look over triangular faces of the 2D mesh
HXTStatus hxtCreateFaceSearchStructure(HXTMesh* mesh, uint32_t **pfaces);
//// creates a mesh with all points of the surface mesh
HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions);
/// Compute sizes at vertices of the mesh from existing edges of the tetrahera
HXTStatus hxtComputeMeshSizeFromTrianglesAndLines (HXTMesh* mesh, HXTDelaunayOptions* delOptions);
HXTStatus hxtComputeMeshSizeFromMesh (HXTMesh* mesh, HXTDelaunayOptions* delOptions);
/// Gives a unique color to each enclosed volume
HXTStatus hxtColorMesh(HXTMesh* mesh, uint16_t *nbColors);
/// Recover the boundary
HXTStatus hxtRecoverBoundary(HXTMesh* mesh);
HXTStatus hxtCreateNodalsizeFromTrianglesAndLines(HXTMesh* mesh, HXTDelaunayOptions* delOptions);
HXTStatus hxtCreateNodalsizeFromMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions);
HXTStatus hxtDestroyNodalsize(HXTDelaunayOptions* delOptions);
/// Add points at tets circumcenter in order to fullfill a mesh size constraint
HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions, HXTMeshSize* meshsize);
......
#include "hxt_mesh3d.h"
#include "hxt_tetrahedra.h"
#include "hxt_tetDelaunay.h"
#include "hxt_tetRepair.h"
#include "hxt_tetUtils.h"
#include "hxt_tetFlag.h"
#include "hxt_tetColor.h"
#include "hxt_tetOpti.h"
HXTStatus hxtTetMesh3d(HXTMesh* mesh,
int nthreads,
int defaulThreads,
int DelaunayThreads,
int optimizationThreads,
int reproducible,
int verbosity,
int displayStat,
int refine, // refine if !=0
int optimize,// optimize quality if !=0
int refine,
int optimize,
double qualityThreshold,
HXTStatus (*bnd_recovery)(HXTMesh* mesh)) {
if(defaulThreads>0) {
omp_set_num_threads(defaulThreads);
}
double t[8]={0};
t[0] = omp_get_wtime();
......@@ -21,7 +30,7 @@ HXTStatus hxtTetMesh3d(HXTMesh* mesh,
hxtBboxInit(&bbox);
HXT_CHECK( hxtBboxAdd(&bbox, mesh->vertices.coord, mesh->vertices.num) );
HXTDelaunayOptions delOptions = {&bbox, NULL, 0.0, 0.0, 0, verbosity, reproducible, nthreads};
HXTDelaunayOptions delOptions = {&bbox, NULL, 0.0, 0.0, 0, verbosity, reproducible, DelaunayThreads};
uint32_t numVerticesConstrained = mesh->vertices.num;
HXT_INFO_COND(verbosity>0, "Creating an empty mesh with %u vertices", numVerticesConstrained);
......@@ -29,22 +38,38 @@ HXTStatus hxtTetMesh3d(HXTMesh* mesh,
t[1] = omp_get_wtime();
uint64_t nbMissingFacets, nbMissingEdges=0;
uint64_t nbMissingTriangles, nbLinesNotInTriangles, nbMissingLines=0;
uint16_t nbColors;
HXT_CHECK( hxtConstrainTriangles(mesh, &nbMissingFacets) );
if(nbMissingFacets==0) // TODO: differentiating missing triangles and missing edges ??
HXT_CHECK( hxtConstrainEdgesNotInTriangles(mesh, &nbMissingEdges) );
uint64_t* tri2TetMap = NULL;
uint64_t* lines2TriMap = NULL;
uint64_t* lines2TetMap = NULL;
HXT_CHECK( hxtAlignedMalloc(&tri2TetMap, mesh->triangles.num*sizeof(uint64_t)) );
HXT_CHECK( hxtAlignedMalloc(&lines2TriMap, mesh->lines.num*sizeof(uint64_t)) );
HXT_CHECK( hxtGetTri2TetMap(mesh, tri2TetMap, &nbMissingTriangles) );
HXT_CHECK( hxtGetLines2TriMap(mesh, lines2TriMap, &nbLinesNotInTriangles) );
if(nbLinesNotInTriangles!=0) {
HXT_CHECK( hxtAlignedMalloc(&lines2TetMap, mesh->lines.num*sizeof(uint64_t)) );
if(nbMissingTriangles==0) {
HXT_CHECK( hxtGetLines2TetMap(mesh, lines2TetMap, &nbMissingLines) );
}
}
t[2] = omp_get_wtime();
if (nbMissingFacets != 0 || nbMissingEdges!=0){
if (nbMissingTriangles!=0 || nbMissingLines!=0){
if(bnd_recovery==NULL)
return HXT_ERROR_MSG(HXT_STATUS_ERROR,
"there are missing features but no boundary recovery function is given");
if(nbMissingFacets)
HXT_INFO("Recovering %lu missing facet(s)", nbMissingFacets);
else if(nbMissingEdges)
HXT_INFO("Recovering %lu missing edge(s)", nbMissingEdges);
if(nbMissingTriangles)
HXT_INFO("Recovering %lu missing facet(s)", nbMissingTriangles);
else if(nbMissingLines)
HXT_INFO("Recovering %lu missing edge(s)", nbMissingLines);
HXT_CHECK(bnd_recovery(mesh));
if(delOptions.numVerticesInMesh < mesh->vertices.num) {
......@@ -54,113 +79,62 @@ HXTStatus hxtTetMesh3d(HXTMesh* mesh,
t[3] = omp_get_wtime();
memset(mesh->tetrahedra.flag, 0, mesh->tetrahedra.num*sizeof(uint16_t));
HXT_CHECK(hxtTetOrientNodes(mesh));
HXT_CHECK(hxtTetAdjacencies(mesh));
HXT_CHECK(hxtAddGhosts(mesh));
// HXT_CHECK( hxtTetVerify(mesh) );
HXT_CHECK( hxtConstrainTriangles(mesh, &nbMissingFacets) );
if(nbMissingFacets!=0)
HXT_CHECK( hxtGetTri2TetMap(mesh, tri2TetMap, &nbMissingTriangles) );
if(nbMissingTriangles!=0)
return HXT_ERROR_MSG( HXT_STATUS_ERROR,
"%d boundary face%s still missing (after recovery step).",
nbMissingFacets, (nbMissingFacets>1)?"s are":" is" );
nbMissingTriangles, (nbMissingTriangles>1)?"s are":" is" );
HXT_CHECK( hxtConstrainEdgesNotInTriangles(mesh, &nbMissingEdges) );
if(nbMissingEdges!=0)
if(nbLinesNotInTriangles!=0)
HXT_CHECK( hxtGetLines2TetMap(mesh, lines2TetMap, &nbMissingLines) );
if(nbMissingLines!=0)
return HXT_ERROR_MSG( HXT_STATUS_ERROR,
"%d constrained edge%s still missing (after recovery step).",
nbMissingEdges, (nbMissingEdges>1)?"s are":" is" );
nbMissingLines, (nbMissingLines>1)?"s are":" is" );
}
HXT_CHECK(hxtColorMesh(mesh, &nbColors));
#ifdef DEBUG
HXT_CHECK( hxtTetVerify(mesh) );
HXT_CHECK( hxtConstrainTriangles(mesh, tri2TetMap) );
if(nbLinesNotInTriangles!=0)
HXT_CHECK( hxtConstrainLinesNotInTriangles(mesh, lines2TetMap, lines2TriMap) );
memset(mesh->tetrahedra.flag, 0, mesh->tetrahedra.num*sizeof(uint16_t));
HXT_CHECK( hxtConstrainTriangles(mesh, &nbMissingFacets) );
if(nbMissingFacets!=0)
return HXT_ERROR_MSG( HXT_STATUS_ERROR,
"%d boundary face%s still missing (after refine).",
nbMissingFacets, (nbMissingFacets>1)?"s are":" is" );
HXT_CHECK( hxtColorMesh(mesh, &nbColors) );
HXT_CHECK( hxtMapColorsToBrep(mesh, nbColors, tri2TetMap) );
HXT_CHECK( hxtConstrainEdgesNotInTriangles(mesh, &nbMissingEdges) );
if(nbMissingEdges!=0)
return HXT_ERROR_MSG( HXT_STATUS_ERROR,
"%d constrained edge%s still missing (after refine).",
nbMissingEdges, (nbMissingEdges>1)?"s are":" is" );
#endif
HXT_CHECK( hxtAlignedFree(&tri2TetMap) );
HXT_CHECK( hxtAlignedFree(&lines2TetMap) );
HXT_CHECK( hxtAlignedFree(&lines2TriMap) );
t[4] = omp_get_wtime();
if(refine){
// HXT_CHECK(hxtComputeMeshSizeFromMesh(mesh, &delOptions));
HXT_CHECK(hxtComputeMeshSizeFromTrianglesAndLines(mesh, &delOptions));
// // triangulate only one color
// if(color>=0) {
// #pragma omp parallel for simd
// for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
// if(mesh->tetrahedra.colors[i]!=theColor){
// markTetAsProcessed(mesh,i) = 1;
// }
// }
// }
HXT_CHECK(hxtCreateNodalsizeFromTrianglesAndLines(mesh, &delOptions));
if(nbColors!=mesh->brep.numVolumes) {
HXT_CHECK( setFlagsToProcessOnlyVolumesInBrep(mesh) );
}
HXTMeshSize *meshSize = NULL;
// HXT_CHECK(hxtMeshSizeCreate (context,&meshSize));
// HXT_CHECK(hxtMeshSizeCompute (meshSize, bbox.min, bbox.max, mySize, NULL));
// printf("time from empty mesh to first insertion: %f second\n", omp_get_wtime() - time);
HXT_CHECK(hxtRefineTetrahedra(mesh, &delOptions, meshSize));
// HXT_CHECK(hxtMeshSizeDelete (&meshSize));
HXT_CHECK(hxtAlignedFree(&delOptions.nodalSizes));
HXT_CHECK( hxtDestroyNodalsize(&delOptions) );
// #endif
}
t[5] = omp_get_wtime();
#ifdef DEBUG
HXT_CHECK( hxtTetVerify(mesh) );
memset(mesh->tetrahedra.flag, 0, mesh->tetrahedra.num*sizeof(uint16_t));
HXT_CHECK( hxtConstrainTriangles(mesh, &nbMissingFacets) );
if(nbMissingFacets!=0)
return HXT_ERROR_MSG( HXT_STATUS_ERROR,
"%d boundary face%s still missing (after refine).",
nbMissingFacets, (nbMissingFacets>1)?"s are":" is" );
HXT_CHECK( hxtConstrainEdgesNotInTriangles(mesh, &nbMissingEdges) );
if(nbMissingEdges!=0)
return HXT_ERROR_MSG( HXT_STATUS_ERROR,
"%d constrained edge%s still missing (after refine).",
nbMissingEdges, (nbMissingEdges>1)?"s are":" is" );
#endif
if(optimize)
HXT_CHECK( hxtOptimizeTetrahedra(mesh, &bbox, delOptions.minSizeEnd, qualityThreshold, numVerticesConstrained) );
HXT_CHECK( hxtOptimizeTetrahedra(mesh, &bbox, optimizationThreads, delOptions.minSizeEnd, qualityThreshold, numVerticesConstrained) );
t[6] = omp_get_wtime();
#ifdef DEBUG
HXT_CHECK( hxtTetVerify(mesh) );
memset(mesh->tetrahedra.flag, 0, mesh->tetrahedra.num*sizeof(uint16_t));
HXT_CHECK( hxtConstrainTriangles(mesh, &nbMissingFacets) );
if(nbMissingFacets!=0)
return HXT_ERROR_MSG( HXT_STATUS_ERROR,
"%d boundary face%s still missing (after refine).",
nbMissingFacets, (nbMissingFacets>1)?"s are":" is" );
HXT_CHECK( hxtConstrainEdgesNotInTriangles(mesh, &nbMissingEdges) );
if(nbMissingEdges!=0)
return HXT_ERROR_MSG( HXT_STATUS_ERROR,
"%d constrained edge%s still missing (after refine).",
nbMissingEdges, (nbMissingEdges>1)?"s are":" is" );
#endif
if(displayStat){
HXT_INFO("\n\t\tFinal tet. mesh contains %lu tetrahedra"
......
......@@ -4,7 +4,9 @@
#include "hxt_mesh.h"
HXTStatus hxtTetMesh3d(HXTMesh* mesh,
int nthreads,
int defaulThreads,
int DelaunayThreads,
int optimizationThreads,
int reproducible,
int verbosity,
int displayStat,
......
......@@ -4,12 +4,13 @@
#include <omp.h>
#else
#include <time.h>