toto

parent 3ed270bc
......@@ -80,6 +80,7 @@ struct contextGeometryOptions {
int copyMeshingMethod, exactExtrusion;
int matchGeomAndMesh;
int hideCompounds, orientedPhysicals, doubleClickedEntityTag;
int reparamOnFaceRobust;
std::string doubleClickedPointCommand, doubleClickedLineCommand;
std::string doubleClickedSurfaceCommand, doubleClickedVolumeCommand;
};
......
......@@ -899,6 +899,10 @@ StringXNumber GeometryOptions_Number[] = {
{ F|O, "PointType" , opt_geometry_point_type , 0. ,
"Display points as solid color dots (0) or 3D spheres (1)" },
{ F|O, "ReparamOnFaceRobust" , opt_geometry_reparam_on_face_robust, 0 ,
"Use projection for reparametrization of a point classified on GEdge on a GFace" },
{ F|O, "ScalingFactor" , opt_geometry_scaling_factor , 1.0 ,
"Global geometry scaling factor" },
{ F|O, "OrientedPhysicals" , opt_geometry_oriented_physicals, 1. ,
......
......@@ -55,6 +55,8 @@ typedef unsigned long intptr_t;
int GmshInitialize(int argc, char **argv)
{
Msg::SetNumThreads(1);
static bool isInitialized = false;
if(isInitialized) return 1;
isInitialized = true;
......
......@@ -775,7 +775,7 @@ void Msg::ProgressMeter(int n, int N, bool log, const char *fmt, ...)
}
#if defined(_OPENMP)
#pragma omp critical
#pragma omp barrier
#endif
}
}
......@@ -1601,4 +1601,4 @@ void MsgProgressStatus::next()
Msg::ProgressMeter(currentI_-1, totalElementToTreat_, true, "%d%% (remaining time ~%g hours)",
currentPercentage, remaining/3600);
}
}
\ No newline at end of file
}
......@@ -18,6 +18,10 @@
#include "StringUtils.h"
#include "Context.h"
#if defined(_OPENMP)
#include <omp.h>
#endif
#if defined(HAVE_ZIPPER)
#include <iostream>
#include <fstream>
......@@ -336,10 +340,14 @@ void CheckResources()
double Cpu()
{
#if defined(_OPENMP)
return omp_get_wtime();
#else
long mem = 0;
double s = 0.;
GetResources(&s, &mem);
return s;
#endif
}
double TotalRam()
......
......@@ -2830,6 +2830,14 @@ double opt_general_display_border_factor(OPT_ARGS_NUM)
return CTX::instance()->displayBorderFactor;
}
double opt_geometry_reparam_on_face_robust(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
CTX::instance()->geom.reparamOnFaceRobust = val;
return CTX::instance()->geom.reparamOnFaceRobust;
}
double opt_general_point_size(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
......
......@@ -376,6 +376,7 @@ double opt_geometry_label_type(OPT_ARGS_NUM);
double opt_geometry_point_size(OPT_ARGS_NUM);
double opt_geometry_point_sel_size(OPT_ARGS_NUM);
double opt_geometry_point_type(OPT_ARGS_NUM);
double opt_geometry_reparam_on_face_robust(OPT_ARGS_NUM);
double opt_geometry_line_width(OPT_ARGS_NUM);
double opt_geometry_line_sel_width(OPT_ARGS_NUM);
double opt_geometry_line_type(OPT_ARGS_NUM);
......
......@@ -121,19 +121,23 @@ SPoint2 OCCEdge::reparamOnFace(const GFace *face, double epar, int dir) const
pnt.Coord(u, v);
// sometimes OCC miserably fails ...
if (0){
if (CTX::instance()->geom.reparamOnFaceRobust){
GPoint p1 = point(epar);
GPoint p2 = face->point(u, v);
const double dx = p1.x()-p2.x();
const double dy = p1.y()-p2.y();
const double dz = p1.z()-p2.z();
if(sqrt(dx * dx + dy * dy + dz * dz) > 1.e-2 * CTX::instance()->lc){
if(sqrt(dx * dx + dy * dy + dz * dz) > CTX::instance()->geom.tolerance){
Msg::Warning("Reparam on face was inaccurate for curve %d on surface %d at point %g",
tag(), face->tag(), epar);
Msg::Warning("On the face %d local (%g %g) global (%g %g %g)",
face->tag(), u, v, p2.x(), p2.y(), p2.z());
Msg::Warning("On the edge %d local (%g) global (%g %g %g)",
tag(), epar, p1.x(), p1.y(), p1.z());
double guess [2] = {u,v};
GPoint pp = face->closestPoint(SPoint3(p1.x(),p1.y(),p1.z()), guess);
u = pp.u();
v = pp.v();
}
}
return SPoint2(u, v);
......
......@@ -91,6 +91,7 @@ class Field {
public:
Field() : update_needed(false) {}
virtual ~Field();
virtual void update() {}
int id;
std::map<std::string, FieldOption *> options;
std::map<std::string, FieldCallback*> callbacks;
......@@ -124,6 +125,7 @@ class FieldManager : public std::map<int, Field*> {
int _boundaryLayer_field;
public:
std::map<std::string, FieldFactory*> map_type_name;
void initialize();
void reset();
Field *get(int id);
Field *newField(int id, std::string type_name);
......
......@@ -37,6 +37,10 @@
#include "meshGRegionRelocateVertex.h"
#include "pointInsertion.h"
#if defined(_OPENMP)
#include <omp.h>
#endif
#if defined(HAVE_OPTHOM)
#include "OptHomRun.h"
#include "OptHomElastic.h"
......@@ -278,6 +282,9 @@ static bool CancelDelaunayHybrid(GModel *m)
static void Mesh0D(GModel *m)
{
m->getFields()->initialize();
for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it){
GVertex *gv = *it;
if(gv->mesh_vertices.empty())
......@@ -299,26 +306,43 @@ static void Mesh0D(GModel *m)
static void Mesh1D(GModel *m)
{
m->getFields()->initialize();
if(TooManyElements(m, 1)) return;
Msg::StatusBar(true, "Meshing 1D...");
double t1 = Cpu();
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
std::vector<GEdge*> temp;
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
(*it)->meshStatistics.status = GEdge::PENDING;
temp.push_back(*it);
}
Msg::ResetProgressMeter();
int nIter = 0, nTot = m->getNumEdges();
while(1){
// meshGEdge mesher;
int nPending = 0;
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
if ((*it)->meshStatistics.status == GEdge::PENDING){
(*it)->mesh(true);
nPending++;
const size_t sss = temp.size();
#if defined(_OPENMP)
#pragma omp parallel for schedule (dynamic)
#endif
for(size_t K = 0 ; K < sss ; K++){
GEdge *ed = temp[K];
if (ed->meshStatistics.status == GEdge::PENDING){
ed->mesh(true);
#if defined(_OPENMP)
#pragma omp critical
#endif
{
nPending++;
}
}
if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 1D...");
// if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 1D...");
}
if(!nPending) break;
if(nIter++ > 10) break;
}
......@@ -394,9 +418,11 @@ static void PrintMesh2dStatistics(GModel *m)
static void Mesh2D(GModel *m)
{
m->getFields()->initialize();
if(TooManyElements(m, 2)) return;
Msg::StatusBar(true, "Meshing 2D...");
double t1 = GetTimeInSeconds();
double t1 = Cpu();
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
(*it)->meshStatistics.status = GFace::PENDING;
......@@ -454,7 +480,7 @@ static void Mesh2D(GModel *m)
nPending++;
}
}
if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D...");
// if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D...");
}
#if defined(_OPENMP)
#pragma omp master
......@@ -487,7 +513,7 @@ static void Mesh2D(GModel *m)
#endif
nPending++;
}
if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D...");
// if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D...");
}
if(!nPending) break;
if(nIter++ > 10) break;
......@@ -496,7 +522,7 @@ static void Mesh2D(GModel *m)
// collapseSmallEdges(*m);
double t2 = GetTimeInSeconds();
double t2 = Cpu();
CTX::instance()->meshTimer[1] = t2 - t1;
Msg::StatusBar(true, "Done meshing 2D (%g s)", CTX::instance()->meshTimer[1]);
......@@ -787,6 +813,8 @@ void TestConformity(GModel *gm)
static void Mesh3D(GModel *m)
{
m->getFields()->initialize();
if(TooManyElements(m, 3)) return;
Msg::StatusBar(true, "Meshing 3D...");
double t1 = Cpu();
......
......@@ -26,10 +26,6 @@
#include "qualityMeasures.h"
#include "meshGRegionLocalMeshMod.h"
#if defined(_HAVE_NUMA)
#include <numa.h>
#endif
//int Tet::in_sphere_counter = 0;
#define SQR(X) (X)*(X)
......@@ -874,26 +870,28 @@ static bool fixDelaunayCavity (Vert *v,
if (_negatives.empty())return false;
// return true;
// unset all tets of the cavity
for (unsigned int i=0; i< cavity.size(); i++)cavity[i]->unset(myThread,K);
for (unsigned int i=0; i<bndK.size(); i++)if(bndK[i].t)bndK[i].t->unset(myThread,K);
// return true;
Msg::Debug("Fixing cavity (%3ld,%3ld) : %ld negatives",
cavity.size(),bndK.size(), _negatives.size());
Tet *containsV = tetContainsV (v,cavity);
if (! containsV) return true;
while (!_negatives.empty()) {
for (unsigned int i=0;i<_negatives.size();i++){
conn &c = bndK[_negatives[i] ];
Tet *toRemove = tetInsideCavityWithFAce (c.f, cavity);
if (toRemove){
std::vector<Tet*>::iterator it = std::find(cavity.begin(), cavity.end(), toRemove);
if (it != cavity.end())
if (it != cavity.end()){
cavity.erase(it);
}
else
Msg::Fatal("Datastructure Broken in %s line %5d",__FILE__,__LINE__);
}
......@@ -902,6 +900,8 @@ static bool fixDelaunayCavity (Vert *v,
buildDelaunayBall (cavity,bndK);
starShapeness (v, bndK, _negatives);
}
for (unsigned int i=0; i< cavity.size(); i++)cavity[i]->set(myThread,K);
for (unsigned int i=0; i<bndK.size(); i++)if(bndK[i].t)bndK[i].t->set(myThread,K);
return false;
}
......@@ -1131,6 +1131,8 @@ void delaunayTrgl (const unsigned int numThreads,
double totCavityGlob=0;
#endif
// printf("%d threads for inserting %d points\n",numThreads,Npts);
// double t1,t2=0,t3=0,t4=0;
// checkLocalDelaunayness(allocator, 0, "initial");
......@@ -1163,7 +1165,7 @@ void delaunayTrgl (const unsigned int numThreads,
std::vector<bool> ok(NPTS_AT_ONCE);
connContainer faceToTet;
std::vector<Tet*> Choice(NPTS_AT_ONCE);
for (unsigned int K=0;K<NPTS_AT_ONCE;K++)Choice[K] = randomTet (myThread, allocator);
for (unsigned int K=0;K<NPTS_AT_ONCE;K++)Choice[K] = randomTet (0, allocator);
invalidCavities [myThread] = 0;
......@@ -1173,12 +1175,7 @@ void delaunayTrgl (const unsigned int numThreads,
for (unsigned int K=0;K<NPTS_AT_ONCE;K++){
locSizeK[K] = assignTo[K+myThread*NPTS_AT_ONCE].size();
locSize += locSizeK[K];
#if defined(_HAVE_NUMA)
allocatedVerts [K] = (Vert*)numa_alloc_local (locSizeK[K]*sizeof(Vert));
#else
// allocatedVerts [K] = (Vert*)calloc (locSizeK[K],sizeof(Vert));
allocatedVerts [K] = new Vert [locSizeK[K]];
#endif
for (unsigned int iP=0 ; iP < locSizeK[K] ; iP++){
allocatedVerts[K][iP] = *(assignTo[K+myThread*NPTS_AT_ONCE][iP]);
if (numThreads!=1) allocatedVerts[K][iP]._thread = myThread;
......@@ -1187,6 +1184,8 @@ void delaunayTrgl (const unsigned int numThreads,
std::vector<Vert*> vToAdd(NPTS_AT_ONCE);
// printf("reaching parallel section\n");
#if defined(_OPENMP)
#pragma omp barrier
#endif
......@@ -1208,7 +1207,7 @@ void delaunayTrgl (const unsigned int numThreads,
vToAdd[K] = iPGlob < locSizeK[K] ? &allocatedVerts[K][iPGlob] : NULL;
if(vToAdd[K]){
// In 3D, insertion of a point may lead to deletion of tets !!
if (!Choice[K]->V[0])Choice[K] = randomTet (myThread, allocator);
if (!Choice[K]->V[0])Choice[K] = randomTet (0, allocator);
// int nbCoucou=0;
while(1){
t[K] = walk ( Choice[K] , vToAdd[K], Npts, totSearch, myThread);
......@@ -1216,7 +1215,7 @@ void delaunayTrgl (const unsigned int numThreads,
// printf("coucou %d\n",nbCoucou++);
// the domain may not be convex. we then start from a random tet and
// walk from there
Choice[K] = randomTet(rand() % numThreads, allocator);
Choice[K] = randomTet(0, allocator);
}
}
}
......@@ -1311,12 +1310,6 @@ void delaunayTrgl (const unsigned int numThreads,
if (invalidCavities[0])Msg::Warning("%d invalid cavities",invalidCavities[0]);
//checkLocalDelaunayness(T,"final");
// printf(" %12.5E %12.5E %12.5E tot %12.5E \n",t2,t3,t4,t2+t3+t4);
#if defined(_VERBOSE)
printf("average searches per point %12.5E\n",totSearchGlob/Npts);
printf("average size for del cavity %12.5E\n",totCavityGlob/Npts);
......@@ -1328,6 +1321,10 @@ void delaunayTrgl (const unsigned int numThreads,
#endif
for (unsigned int myThread=0; myThread < numThreads;myThread++)
for (unsigned int i=0;i<allocator.size(myThread);i++)allocator(myThread,i)->setAllDeleted();
}
static void initialCube (std::vector<Vert*> &v,
......
......@@ -16,7 +16,7 @@
#endif
#ifndef MAX_NUM_THREADS_
#define MAX_NUM_THREADS_ 1
#define MAX_NUM_THREADS_ 8
#endif
typedef unsigned char CHECKTYPE;
......@@ -360,7 +360,10 @@ public:
class tetContainer {
std::vector<aBunchOfStuff<Tet> *> _perThread;
public:
unsigned int size(int thread) const { return _perThread[thread]->size(); }
unsigned int size(int thread) const {
if (_perThread.size() <= thread)return 0;
return _perThread[thread]->size();
}
inline Tet * operator () (int thread, int j) const
{
return (*_perThread[thread])(j);
......@@ -368,7 +371,7 @@ class tetContainer {
tetContainer (int nbThreads, int preallocSizePerThread)
{
// FIXME !!!
if (nbThreads != 1) throw;
// if (nbThreads != 1) throw;
_perThread.resize(nbThreads);
#if defined(_OPENMP)
#pragma omp parallel num_threads(nbThreads)
......
......@@ -199,30 +199,63 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp)
// exit(1);
}
inline bool distributeEdgeThroughThreads ( int nbThreads , int myThread ,
const Edge &ed){
// const size_t h = ((size_t)ed.first->getNum());
const size_t h = ((size_t)ed.first) >> 3;
// printf("%lu %d %d %d\n",h,nbThreads,myThread,h % nbThreads);
return h % nbThreads == myThread;
}
void saturateEdges(edgeContainer &ec,
tetContainer &T,
int nbThreads,
std::vector<Vert*> &S)
{
// _C1 = 0;
// _C2 = 0;
std::stack<IPT> temp;
// FIXME
const int N = T.size(0);
for (int i=0;i<N;i++){
Tet *t = T(0,i);
if (t->V[0] && t->_modified){
t->_modified = false;
for (int iEdge=0;iEdge<6;iEdge++){
Edge ed = t->getEdge(iEdge);
bool isNew = ec.addNewEdge(ed);
if (isNew){
// saturateEdge_A_LA_FACON_PaulLouis(ed, S);
saturateEdge (ed, S, temp);
#pragma omp parallel
{
std::vector<Vert*> Sloc;
std::stack<IPT> temp;
#ifdef _OPENMP
int myThread = omp_get_thread_num();
// nbThreads = omp_get_num_threads();
#else
int myThread = 0;
#endif
for (int iThread = 0; iThread<nbThreads;iThread++){
const int N = T.size(iThread);
for (int i=0;i<N;i++){
Tet *t = T(iThread,i);
if (t->V[0] && t->_modified){
// t->_modified = false;
for (int iEdge=0;iEdge<6;iEdge++){
Edge ed = t->getEdge(iEdge);
if (distributeEdgeThroughThreads (nbThreads,myThread, ed)){
bool isNew = ec.addNewEdge(ed);
if (isNew){
saturateEdge (ed, Sloc, temp);
}
}
}
}
}
}
#pragma omp critical
S.insert (S.end(),Sloc.begin(), Sloc.end());
#pragma omp barrier
}
for (int iThread = 0; iThread<nbThreads;iThread++){
const int N = T.size(iThread);
for (int i=0;i<N;i++){
Tet *t = T(iThread,i);
if (t->V[0] && t->_modified){
t->_modified = false;
}
}
}
// printf("timings %12.5E %12.5E\n",_C1,_C2);
}
......@@ -451,6 +484,48 @@ void computeMeshSizes (GRegion *gr, std::map<MVertex*, double> &s){
}
}
void parallelDelaunay (int NT, std::vector<Vert*> &S,
tetContainer &allocator,
int iter,
bool explicitFiltering,
std::vector<int> &indices,
edgeContainer *embeddedEdges){
// Msg::Info ("Parallel Delaunay with %d threads",NT);
int N = S.size();
NT = std::min(NT, MAX_NUM_THREADS_);
std::vector<Vert*> assignTo0[1];
std::vector<std::vector<Vert*> > assignTo (NT);
for (int i=1;i<indices.size();i++){
int start = indices[i-1];
int end = indices[i];
int sizePerBlock = (NT*((end-start) / NT))/NT;
int currentBlock = 0;
int localCounter = 0;
// printf("sizePerBlock[%d] = %d (%d,%d)\n",i,sizePerBlock, start, end);
if (i < 2){
for (int jPt=start;jPt<end;jPt++){
assignTo0[0].push_back(S[jPt]);
}
}
else {
for (int jPt=start;jPt<end;jPt++){
if (localCounter++ >= sizePerBlock && currentBlock != NT-1){
localCounter = 0;
currentBlock++;
}
assignTo[currentBlock].push_back(S[jPt]);
}
}
}
// delaunayTrgl (numThreads,1,add.size(), &add,allocator,embeddedEdges.empty() ? NULL : &embeddedEdges);
delaunayTrgl(1,1, assignTo0[0].size(),assignTo0,allocator,embeddedEdges);
delaunayTrgl(NT,1, N,&assignTo[0], allocator,embeddedEdges);
}
extern double tetQuality (Tet* t, double *volume);
void edgeBasedRefinement(const int numThreads,
......@@ -557,7 +632,7 @@ void edgeBasedRefinement(const int numThreads,
while(1){
std::vector<Vert*> add;
double t1 = Cpu();
// C_COUNT = 0;
// PARALLEL DONE
saturateEdges (ec, allocator, numThreads, add);
// printf("%d points added",add.size());
double t2 = Cpu();
......@@ -572,17 +647,23 @@ void edgeBasedRefinement(const int numThreads,
SortHilbert(add, indices);
// printf("%d indices\n",indices.size());
double t4 = Cpu();
delaunayTrgl (1,1,add.size(), &add,allocator,embeddedEdges.empty() ? NULL : &embeddedEdges);
parallelDelaunay ( numThreads, add, allocator, iter, false, indices,embeddedEdges.empty() ? NULL : &embeddedEdges);
// delaunayTrgl (2 ,1,add.size(), &add,allocator,embeddedEdges.empty() ? NULL : &embeddedEdges);
double t5 = Cpu();
add_all.insert (add_all.end(), add.begin(), add.end());
Msg::Info("IT %3d %8d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",
size_t sss = 0;
for (int myThread=0; myThread < numThreads; myThread++)sss+= allocator.size(myThread);
Msg::Info("IT %3d %8d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",
iter,add.size(),
(t2-t1),
(t3-t2),
(t4-t3),
(t5-t4),
(t5-__t__),
allocator.size(0));
sss);
iter++;
}
}
......@@ -618,26 +699,28 @@ void edgeBasedRefinement(const int numThreads,
int cat [10] = {0,0,0,0,0,0,0,0,0,0};
double MIN = 1.0;
for (unsigned int i=0; i< allocator.size(0);i++){
Tet *tt = allocator (0,i);
MVertex *mvs[4];
if (tt->V[0]){
double vol;
double q = tetQuality (tt,&vol);
cat [(int)(q*10)] ++;
MIN = std::min(MIN,q);
for (int j=0;j<4;j++){
Vert *v = tt->V[j];
std::map<Vert*,MVertex*>::iterator it = _ma.find(v);
if (it == _ma.end()){
MVertex *mv = new MVertex (v->x(),v->y(),v->z(),gr);
gr->mesh_vertices.push_back(mv);
_ma[v] = mv;
mvs[j] = mv;
for (int myThread=0; myThread < numThreads; myThread++) {
for (unsigned int i=0; i< allocator.size(myThread);i++){
Tet *tt = allocator (myThread,i);
MVertex *mvs[4];
if (tt->V[0]){
double vol;
double q = tetQuality (tt,&vol);
cat [(int)(q*10)] ++;
MIN = std::min(MIN,q);
for (int j=0;j<4;j++){
Vert *v = tt->V[j];
std::map<Vert*,MVertex*>::iterator it = _ma.find(v);
if (it == _ma.end()){
MVertex *mv = new MVertex (v->x(),v->y(),v->z(),gr);
gr->mesh_vertices.push_back(mv);
_ma[v] = mv;
mvs[j] = mv;
}
else mvs[j] = it->second;
}
else mvs[j] = it->second;
gr->tetrahedra.push_back(new MTetrahedron(mvs[0],mvs[1],mvs[2],mvs[3]));
}
gr->tetrahedra.push_back(new MTetrahedron(mvs[0],mvs[1],mvs[2],mvs[3]));
}
}
......
......@@ -508,8 +508,6 @@ static void addBoundaryLayerPoints(GEdge *ge, double &t_begin, double &t_end,
void meshGEdge::operator() (GEdge *ge)
{
ge->model()->setCurrentMeshEntity(ge);
// if(ge->geomType() == GEntity::DiscreteCurve) return;
if(ge->geomType() == GEntity::BoundaryLayerCurve) return;
if(ge->meshAttributes.method == MESH_NONE) return;
......@@ -534,7 +532,19 @@ void meshGEdge::operator() (GEdge *ge)
return;
}
if(ge->model()->getNumEdges() > 1000){
if(ge->model()->getNumEdges() > 100000){
if (ge->tag() % 100000 == 1){
Msg::Info("Meshing curve %d/%d (%s)", ge->tag(), ge->model()->getNumEdges(),
ge->getTypeString().c_str());
}
}
else if(ge->model()->getNumEdges() > 10000){
if (ge->tag() % 10000 == 1){
Msg::Info("Meshing curve %d/%d (%s)", ge->tag(), ge->model()->getNumEdges(),
ge->getTypeString().c_str());
}
}
else if(ge->model()->getNumEdges() > 1000){
if (ge->tag() % 1000 == 1){
Msg::Info("Meshing curve %d/%d (%s)", ge->tag(), ge->model()->getNumEdges(),
ge->getTypeString().c_str());
......
......@@ -2415,7 +2415,7 @@ void deMeshGFace::operator() (GFace *gf)
}
// for debugging, change value from -1 to -100;
int debugSurface = -1; //-100;
int debugSurface = -100; //-100;
void meshGFace::operator() (GFace *gf, bool print)
{
......
......@@ -63,6 +63,7 @@ static void computeMeshMetricsForBamg(GFace *gf, int numV,
J.mult(M,W);
W.mult(JT,R);
bamg::Metric M1(R(0,0),R(1,0),R(1,1));
// printf("%12.5E %12.5E %12.5E vs %12.5E %12.5E %12.5E \n",m(0,0),m(1,0),m(1,1),M1.a11,M1.a21,M1.a22);
mm11[i] = M1.a11;
mm12[i] = M1.a21;
mm22[i] = M1.a22;
......
......@@ -746,13 +746,13 @@ static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> &regions)
}
else{
void edgeBasedRefinement(const int numThreads,
const int nptsatonce,
GRegion *gr);
const int nptsatonce,
GRegion *gr);
// just to remove tets that are not to be meshed
insertVerticesInRegion(gr, 0);
for(unsigned int i = 0; i < regions.size(); i++){
Msg::Info("Refining volume %d",regions[i]->tag());
edgeBasedRefinement(1, 1, regions[i]);
Msg::Info("Refining volume %d with %d threads",regions[i]->tag(),Msg::GetMaxThreads());
edgeBasedRefinement(Msg::GetMaxThreads(), 1, regions[i]);
}
// RelocateVertices (regions,-1);
}
......
......@@ -825,9 +825,9 @@ bool tetgenmesh::reconstructmesh(void *p)
MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])];
MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])];
MVertex *v3 = pointmark(p[2]) >= _vertices.size() ? _extras[pointmark(p[2])] : _vertices[pointmark(p[2])];
if (pointmark(p[0]) >= _vertices.size())printf("F %d %d\n",v1->getIndex(),pointmark(p[0]));
if (pointmark(p[1]) >= _vertices.size())printf("F %d %d\n",v2->getIndex(),pointmark(p[1]));
if (pointmark(p[2]) >= _vertices.size())printf("F %d %d\n",v3->getIndex(),pointmark(p[2]));
// if (pointmark(p[0]) >= _vertices.size())printf("F %d %d\n",v1->getIndex(),pointmark(p[0]));
// if (pointmark(p[1]) >= _vertices.size())printf("F %d %d\n",v2->getIndex(),pointmark(p[1]));
// if (pointmark(p[2]) >= _vertices.size())printf("F %d %d\n",v3->getIndex(),pointmark(p[2]));
// MVertex *v1 = _vertices[pointmark(p[0])];
// MVertex *v2 = _vertices[pointmark(p[1])];
// MVertex *v3 = _vertices[pointmark(p[2])];
......
......@@ -23,6 +23,43 @@
int MTet4::radiusNorm = 2;