Merge branch 'master' into master

This commit is contained in:
Tim Vincent 2019-11-11 15:19:55 -05:00 committed by GitHub
commit f6c353c359
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
96 changed files with 23088 additions and 3180 deletions

View File

@ -56,6 +56,9 @@ build_win:
-G Ninja
-DCMAKE_INSTALL_PREFIX=%INSTALL_DIR%
-DUSE_OCC=ON
-DOCC_LIBRARY=C:/install_opencascade_7.4.0_static/win64/vc14/lib/TKernel.lib
-DOCC_INCLUDE_DIR=C:/install_opencascade_7.4.0_static/inc
-DOCC_LINK_FREETYPE=ON
-DUSE_CCACHE=ON
-DENABLE_UNIT_TESTS=ON
-DCMAKE_BUILD_TYPE=Release
@ -65,6 +68,7 @@ test_win:
<<: *win
stage: test
script:
- cd tests\pytest
- cd %NETGEN_BUILD_DIR%\netgen
- ctest -C Release -V --output-on-failure
- cd ..
@ -238,6 +242,10 @@ build_mac:
-DENABLE_UNIT_TESTS=ON
-DCMAKE_OSX_DEPLOYMENT_TARGET=10.12
-DCMAKE_OSX_SYSROOT=/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk
-DUSE_OCC=ON
-DOCC_LIBRARY=/usr/local/opt/opencascade-7.4.0/lib/libTKernel.a
-DOCC_INCLUDE_DIR=/usr/local/opt/opencascade-7.4.0/include/opencascade
-DOCC_LINK_FREETYPE=ON
- make -j5 install
test_mac:

View File

@ -14,6 +14,7 @@ get_filename_component(NETGEN_RESOURCE_DIR "${NETGEN_CMAKE_DIR}/@NETGEN_RESOURCE
set(NETGEN_SOURCE_DIR "@PROJECT_SOURCE_DIR@")
set(NETGEN_CHECK_RANGE "@CHECK_RANGE@")
set(NETGEN_INCLUDE_DIRS "${NETGEN_INCLUDE_DIR}/include;${NETGEN_INCLUDE_DIR}")
set(NETGEN_CMAKE_THREAD_LIBS_INIT "@CMAKE_THREAD_LIBS_INIT@")
set(NETGEN_FFMPEG_LIBRARIES "@FFMPEG_LIBRARIES@")

View File

@ -80,6 +80,14 @@ set(OCC_LIBRARY_NAMES
TKXSBase
)
if(OCC_LINK_FREETYPE)
set(OCC_LIBRARY_NAMES ${OCC_LIBRARY_NAMES} freetype)
endif(OCC_LINK_FREETYPE)
if(OCC_VERSION_STRING VERSION_GREATER_EQUAL "7.3.0")
set(OCC_LIBRARY_NAMES ${OCC_LIBRARY_NAMES} TKVCAF)
endif()
foreach( libname ${OCC_LIBRARY_NAMES} )
find_library( ${libname} ${libname} ${OCC_LIBRARY_DIR} )
set(OCC_LIBRARIES ${OCC_LIBRARIES} ${${libname}})

@ -1 +1 @@
Subproject commit 4ffa6cd2d4a277b1cc8bc10d5c73cc0ee8edfaeb
Subproject commit 7ec2ddfc95f65d1e986d359466a6c254aa514ef3

View File

@ -32,6 +32,7 @@ if(WIN32)
endif(WIN32)
target_compile_definitions(ngcore PUBLIC $<$<CONFIG:DEBUG>:NETGEN_ENABLE_CHECK_RANGE>)
target_include_directories(ngcore INTERFACE $<INSTALL_INTERFACE:${NG_INSTALL_DIR_INCLUDE}> $<INSTALL_INTERFACE:${NG_INSTALL_DIR_INCLUDE}/include>)
if(CHECK_RANGE)
target_compile_definitions(ngcore PUBLIC NETGEN_ENABLE_CHECK_RANGE)

View File

@ -14,6 +14,7 @@
#include <memory>
#include <cxxabi.h>
#include <signal.h>
#include <vector>
namespace ngcore
{
@ -104,11 +105,11 @@ namespace ngcore
if(!funcname.empty())
{
std::array<char, 256> buffer;
std::vector<char> buffer(10240);
int status;
size_t size = buffer.size();
abi::__cxa_demangle(funcname.c_str(), buffer.data(), &size, &status);
out << "in " << yellow << buffer.data() << reset_shell << '\n';
abi::__cxa_demangle(funcname.c_str(), &buffer[0], &size, &status);
out << "in " << yellow << &buffer[0] << reset_shell << '\n';
std::string nm_command = "nm " + libname + " | grep " + funcname + " | cut -f 1 -d ' '";
std::string output;
@ -145,12 +146,12 @@ namespace ngcore
{
std::cerr << "Collecting backtrace..." << std::endl;
std::stringstream result;
void *bt[1024];
void *bt[100];
int bt_size;
char **bt_syms;
int i;
bt_size = backtrace(bt, 1024);
bt_size = backtrace(bt, 100);
bt_syms = backtrace_symbols(bt, bt_size);
Dl_info info;
for (i = 1; i < bt_size-1; i++)
@ -167,6 +168,11 @@ namespace ngcore
static void ngcore_signal_handler(int sig)
{
static bool first_call = true;
if(!first_call)
exit(1); // avoid endless recursions if signals are caused by this handler
first_call = false;
switch(sig)
{
case SIGABRT:

View File

@ -827,6 +827,12 @@ namespace ngcore
return ost;
}
template <typename TI>
NETGEN_INLINE size_t HashValue (const INT<3,TI> ind)
{
INT<3,size_t> lind = ind;
return 113*lind[0]+59*lind[1]+lind[2];
}
template <typename TI>
NETGEN_INLINE size_t HashValue (const INT<2,TI> ind)

View File

@ -15,7 +15,7 @@ namespace ngcore
{
std::ostream* testout = new std::ostream(nullptr); // NOLINT
level::level_enum Logger::global_level;
level::level_enum Logger::global_level = level::warn;
void Logger::log(level::level_enum level, std::string && s)
{

View File

@ -74,6 +74,22 @@ inline void operator delete[]( void* ptr, std::align_val_t al ) noexcept
delete[] (char*)ptr;
}
inline void operator delete ( void* ptr, std::size_t sz, std::align_val_t al ) noexcept
{
if (int(al) > __STDCPP_DEFAULT_NEW_ALIGNMENT__)
_mm_free(ptr);
else
delete (char*)ptr;
}
inline void operator delete[]( void* ptr, std::size_t sz, std::align_val_t al ) noexcept
{
if (int(al) > __STDCPP_DEFAULT_NEW_ALIGNMENT__)
_mm_free(ptr);
else
delete[] (char*)ptr;
}
#endif // __MAC_OS_X_VERSION_MIN_REQUIRED
#endif // __MAC_OS_X_VERSION_MIN_REQUIRED < 101300

View File

@ -35,11 +35,11 @@ namespace ngcore
int TaskManager :: num_threads = 1;
#ifndef __clang__
// #ifndef __clang__
thread_local int TaskManager :: thread_id = 0;
#else
__thread int TaskManager :: thread_id;
#endif
// #else
// __thread int TaskManager :: thread_id;
// #endif
const function<void(TaskInfo&)> * TaskManager::func;
const function<void()> * TaskManager::startup_function = nullptr;

View File

@ -15,6 +15,7 @@
#include "array.hpp"
#include "paje_trace.hpp"
#include "profiler.hpp"
namespace ngcore
{
@ -71,11 +72,11 @@ namespace ngcore
#ifndef __clang__
// #ifndef __clang__
static thread_local int thread_id;
#else
static __thread int thread_id;
#endif
// #else
// static __thread int thread_id;
// #endif
NGCORE_API static bool use_paje_trace;
public:
@ -1059,6 +1060,7 @@ public:
template <typename Tmask>
int ComputeColoring( FlatArray<int> colors, size_t ndofs, Tmask const & getDofs)
{
static Timer timer("ComputeColoring - "+Demangle(typeid(Tmask).name())); RegionTimer rt(timer);
static_assert(sizeof(unsigned int)==4, "Adapt type of mask array");
auto n = colors.Size();

View File

@ -15,10 +15,22 @@ namespace ngcore
// windows does demangling in typeid(T).name()
NGCORE_API std::string Demangle(const char* typeinfo) { return typeinfo; }
#else
NGCORE_API std::string Demangle(const char* typeinfo) { int status; return abi::__cxa_demangle(typeinfo,
nullptr,
nullptr,
&status); }
NGCORE_API std::string Demangle(const char* typeinfo)
{
int status=0;
try
{
char *s = abi::__cxa_demangle(typeinfo, nullptr, nullptr, &status);
std::string result{s};
free(s);
return result;
}
catch( const std::exception & e )
{
GetLogger("utils")->warn("{}:{} cannot demangle {}, status: {}, error:{}", __FILE__, __LINE__, typeinfo, status, e.what());
}
return typeinfo;
}
#endif
double seconds_per_tick = [] () noexcept

View File

@ -72,6 +72,92 @@ namespace netgen
Clean();
}
PointGeomInfo CSGeometry :: ProjectPoint(int surfind, Point<3> & p) const
{
Point<3> hp = p;
GetSurface(surfind)->Project (hp);
p = hp;
return PointGeomInfo();
}
bool CSGeometry :: ProjectPointGI(int surfind, Point<3> & p, PointGeomInfo & gi) const
{
GetSurface(surfind)->Project (p);
return true;
}
void CSGeometry :: ProjectPointEdge(int surfind, INDEX surfind2,
Point<3> & p, EdgePointGeomInfo* /*unused*/) const
{
Point<3> hp = p;
ProjectToEdge (GetSurface(surfind),
GetSurface(surfind2), hp);
p = hp;
}
Vec<3> CSGeometry :: GetNormal(int surfind, const Point<3> & p,
const PointGeomInfo* /*unused*/) const
{
Vec<3> hn;
GetSurface(surfind)->CalcGradient(p, hn);
hn.Normalize();
return hn;
}
void CSGeometry ::
PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi != -1)
{
GetSurface (surfi) -> Project (hnewp);
newgi.trignum = 1;
}
newp = hnewp;
}
void CSGeometry :: PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
Point<3> hnewp = p1+secpoint*(p2-p1);
//(*testout) << "hnewp " << hnewp << " s1 " << surfi1 << " s2 " << surfi2 << endl;
if (surfi1 != -1 && surfi2 != -1 && surfi1 != surfi2)
{
netgen::ProjectToEdge (GetSurface(surfi1),
GetSurface(surfi2),
hnewp);
// (*testout) << "Pointbetween, newp = " << hnewp << endl
// << ", err = " << sqrt (sqr (hnewp(0))+ sqr(hnewp(1)) + sqr (hnewp(2))) - 1 << endl;
newgi.edgenr = 1;
//(*testout) << "hnewp (a1) " << hnewp << endl;
}
else if (surfi1 != -1)
{
GetSurface (surfi1) -> Project (hnewp);
//(*testout) << "hnewp (a2) " << hnewp << endl;
}
newp = hnewp;
};
Vec<3> CSGeometry :: GetTangent(const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
Vec<3> n1 = GetSurface (surfi1)->GetNormalVector (p);
Vec<3> n2 = GetSurface (surfi2)->GetNormalVector (p);
Vec<3> tau = Cross (n1, n2).Normalize();
return tau;
}
void CSGeometry :: Clean ()
{
@ -137,15 +223,6 @@ namespace netgen
return CSGGenerateMesh (*this, mesh, mparam);
}
const Refinement & CSGeometry :: GetRefinement () const
{
// cout << "get CSGeometry - Refinement" << endl;
// should become class variables
RefinementSurfaces * ref = new RefinementSurfaces(*this);
ref -> Set2dOptimizer(new MeshOptimize2dSurfaces(*this));
return *ref;
}
class WritePrimitivesIt : public SolidIterator
{
ostream & ost;

View File

@ -188,6 +188,27 @@ namespace netgen
virtual void SaveToMeshFile (ostream & ost) const override;
PointGeomInfo ProjectPoint(INDEX surfind, Point<3> & p) const override;
bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
void ProjectPointEdge(INDEX surfind, INDEX surfind2, Point<3> & p,
EdgePointGeomInfo* gi = nullptr) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const override;
void PointBetween(const Point<3> & p1, const Point<3> & p2,
double secpoint, int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const override;
int GetChangeVal() { return changeval; }
void Change() { changeval++; }
@ -348,8 +369,6 @@ namespace netgen
virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam) override;
virtual const Refinement & GetRefinement () const override;
void AddSplineSurface (shared_ptr<SplineSurface> ss) { spline_surfaces.Append(ss); }
};

View File

@ -422,7 +422,7 @@ namespace netgen
geom.GetSurface((mesh.GetFaceDescriptor(k).SurfNr()));
Meshing2Surfaces meshing(*surf, mparam, geom.BoundingBox());
Meshing2Surfaces meshing(geom, *surf, mparam, geom.BoundingBox());
meshing.SetStartTime (starttime);
double eps = 1e-8 * geom.MaxSize();
@ -523,48 +523,48 @@ namespace netgen
if (multithread.terminate) return;
{
MeshOptimize2dSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2));
meshopt.EdgeSwapping (i > mparam.optsteps2d/2);
}
if (multithread.terminate) return;
{
// mesh.CalcSurfacesOfNode();
MeshOptimize2dSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mesh, mparam);
meshopt.ImproveMesh(mparam);
}
{
MeshOptimize2dSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.CombineImprove (mesh);
meshopt.CombineImprove();
// mesh.CalcSurfacesOfNode();
}
if (multithread.terminate) return;
{
MeshOptimize2dSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mesh, mparam);
meshopt.ImproveMesh(mparam);
}
}
}

View File

@ -14,13 +14,14 @@ Meshing2Surfaces :: Meshing2Surfaces (const Surface & asurface)
;
}
*/
Meshing2Surfaces :: Meshing2Surfaces (const Surface & asurf,
Meshing2Surfaces :: Meshing2Surfaces (const CSGeometry& geo,
const Surface & asurf,
const MeshingParameters & mp,
const Box<3> & abb)
: Meshing2(mp, abb), surface(asurf), mparam (mp)
{
: Meshing2(geo, mp, abb), surface(asurf), mparam (mp)
{
;
}
}
void Meshing2Surfaces :: DefineTransformation (const Point<3> & p1, const Point<3> & p2,
@ -59,147 +60,4 @@ double Meshing2Surfaces :: CalcLocalH (const Point<3> & p, double gh) const
return loch;
*/
}
MeshOptimize2dSurfaces :: MeshOptimize2dSurfaces (const CSGeometry & ageometry)
: MeshOptimize2d(), geometry(ageometry)
{
;
}
void MeshOptimize2dSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
Point<3> hp = p;
geometry.GetSurface(surfind)->Project (hp);
p = hp;
}
void MeshOptimize2dSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2,
Point<3> & p) const
{
Point<3> hp = p;
ProjectToEdge ( geometry.GetSurface(surfind),
geometry.GetSurface(surfind2), hp);
p = hp;
}
void MeshOptimize2dSurfaces ::
GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const
{
Vec<3> hn = n;
geometry.GetSurface(surfind)->CalcGradient (p, hn);
hn.Normalize();
n = hn;
/*
if (geometry.GetSurface(surfind)->Inverse())
n *= -1;
*/
}
RefinementSurfaces :: RefinementSurfaces (const CSGeometry & ageometry)
: Refinement(), geometry(ageometry)
{
if(geometry.GetNSurf() == 0)
*testout << endl
<< "WARNING: Initializing 2D refinement with 0-surface geometry" << endl
<< "==========================================================" << endl
<< endl << endl;
}
RefinementSurfaces :: ~RefinementSurfaces ()
{
;
}
void RefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi != -1)
{
geometry.GetSurface (surfi) -> Project (hnewp);
newgi.trignum = 1;
}
newp = hnewp;
}
void RefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
Point<3> hnewp = p1+secpoint*(p2-p1);
//(*testout) << "hnewp " << hnewp << " s1 " << surfi1 << " s2 " << surfi2 << endl;
if (surfi1 != -1 && surfi2 != -1 && surfi1 != surfi2)
{
netgen::ProjectToEdge (geometry.GetSurface(surfi1),
geometry.GetSurface(surfi2),
hnewp);
// (*testout) << "Pointbetween, newp = " << hnewp << endl
// << ", err = " << sqrt (sqr (hnewp(0))+ sqr(hnewp(1)) + sqr (hnewp(2))) - 1 << endl;
newgi.edgenr = 1;
//(*testout) << "hnewp (a1) " << hnewp << endl;
}
else if (surfi1 != -1)
{
geometry.GetSurface (surfi1) -> Project (hnewp);
//(*testout) << "hnewp (a2) " << hnewp << endl;
}
newp = hnewp;
};
Vec<3> RefinementSurfaces :: GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
Vec<3> n1 = geometry.GetSurface (surfi1)->GetNormalVector (p);
Vec<3> n2 = geometry.GetSurface (surfi2)->GetNormalVector (p);
Vec<3> tau = Cross (n1, n2).Normalize();
return tau;
}
Vec<3> RefinementSurfaces :: GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const
{
return geometry.GetSurface (surfi1)->GetNormalVector (p);
}
void RefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi) const
{
if (surfi != -1)
geometry.GetSurface (surfi) -> Project (p);
};
void RefinementSurfaces :: ProjectToEdge (Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const
{
netgen::ProjectToEdge (geometry.GetSurface(surfi1),
geometry.GetSurface(surfi2),
p);
}
}

View File

@ -16,7 +16,9 @@ namespace netgen
///
// Meshing2Surfaces (const Surface & asurf);
///
Meshing2Surfaces (const Surface & asurf, const MeshingParameters & mp,
Meshing2Surfaces (const CSGeometry& geo,
const Surface & asurf,
const MeshingParameters & mp,
const Box<3> & aboundingbox);
protected:
@ -38,62 +40,6 @@ namespace netgen
///
double CalcLocalH(const Point<3> & p, double gh) const override;
};
///
class MeshOptimize2dSurfaces : public MeshOptimize2d
{
///
const CSGeometry & geometry;
public:
///
MeshOptimize2dSurfaces (const CSGeometry & ageometry);
///
virtual void ProjectPoint (INDEX surfind, Point<3> & p) const override;
///
virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const override;
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const override;
};
class RefinementSurfaces : public Refinement
{
const CSGeometry & geometry;
public:
RefinementSurfaces (const CSGeometry & ageometry);
virtual ~RefinementSurfaces ();
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
virtual Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const override;
virtual Vec<3> GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi) const override;
virtual void ProjectToEdge (Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const override;
};
}
#endif

View File

@ -1390,7 +1390,6 @@ inline size_t HashValue (INDEX_2 i2, size_t size) { return (113*size_t(i2[0])+si
ClosedHashTable (ClosedHashTable && ht2) = default;
// who needs that ?
ClosedHashTable (NgFlatArray<T_HASH> _hash, NgFlatArray<T> _cont)
: size(_hash.Size()), used(0), hash(_hash.Size(), _hash.Addr(0)), cont(_cont.Size(), _cont.Addr(0))
{

View File

@ -96,12 +96,6 @@ void ParallelFor( int first, int next, const TFunc & f )
template<typename T>
inline atomic<T> & AsAtomic (T & d)
{
return reinterpret_cast<atomic<T>&> (d);
}
typedef void (*TaskManager)(std::function<void(int,int)>);
typedef void (*Tracer)(string, bool); // false .. start, true .. stop

View File

@ -1,5 +1,5 @@
add_definitions(-DNGLIB_EXPORTS)
add_library(geom2d ${NG_LIB_TYPE} genmesh2d.cpp geom2dmesh.cpp geometry2d.cpp python_geom2d.cpp )
add_library(geom2d ${NG_LIB_TYPE} genmesh2d.cpp geometry2d.cpp python_geom2d.cpp )
if(APPLE)
set_target_properties( geom2d PROPERTIES SUFFIX ".so")
endif(APPLE)
@ -18,7 +18,7 @@ if(USE_GUI)
endif(USE_GUI)
install(FILES
geom2dmesh.hpp geometry2d.hpp spline2d.hpp
geometry2d.hpp spline2d.hpp
vsgeom2d.hpp
DESTINATION ${NG_INSTALL_DIR_INCLUDE}/geom2d COMPONENT netgen_devel
)

View File

@ -565,7 +565,7 @@ namespace netgen
mp.quad = hquad || geometry.GetDomainQuadMeshing (domnr);
Meshing2 meshing (mp, Box<3> (pmin, pmax));
Meshing2 meshing (geometry, mp, Box<3> (pmin, pmax));
NgArray<int, PointIndex::BASE> compress(bnp);
compress = -1;

View File

@ -1,119 +0,0 @@
#include <meshing.hpp>
#include <geometry2d.hpp>
namespace netgen
{
Refinement2d :: Refinement2d (const SplineGeometry2d & ageometry)
: Refinement(), geometry(ageometry)
{
;
}
Refinement2d :: ~Refinement2d ()
{
;
}
void Refinement2d ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
newgi.trignum = 1;
}
void Refinement2d ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
Point<2> p2d;
double newdist;
auto spline = geometry.GetSplines().Get(ap1.edgenr);
if( (ap1.dist == 0.0) && (ap2.dist == 0.0) )
{
// used for manually generated meshes
const SplineSeg3<2> * ss3;
const LineSeg<2> * ls;
auto ext = dynamic_cast<const SplineSegExt *>(spline);
if(ext)
{
ss3 = dynamic_cast<const SplineSeg3<2> *>(ext->seg);
ls = dynamic_cast<const LineSeg<2> *>(ext->seg);
}
else
{
ss3 = dynamic_cast<const SplineSeg3<2> *>(spline);
ls = dynamic_cast<const LineSeg<2> *>(spline);
}
Point<2> p12d(p1(0),p1(1)), p22d(p2(0),p2(1));
Point<2> p1_proj(0.0,0.0), p2_proj(0.0,0.0);
double t1_proj = 0.0;
double t2_proj = 0.0;
if(ss3)
{
ss3->Project(p12d,p1_proj,t1_proj);
ss3->Project(p22d,p2_proj,t2_proj);
}
else if(ls)
{
ls->Project(p12d,p1_proj,t1_proj);
ls->Project(p22d,p2_proj,t2_proj);
}
p2d = spline->GetPoint (((1-secpoint)*t1_proj+secpoint*t2_proj));
newdist = (1-secpoint)*t1_proj+secpoint*t2_proj;
}
else
{
p2d = spline->GetPoint (((1-secpoint)*ap1.dist+secpoint*ap2.dist));
newdist = (1-secpoint)*ap1.dist+secpoint*ap2.dist;
}
// (*testout) << "refine 2d line, ap1.dist, ap2.dist = " << ap1.dist << ", " << ap2.dist << endl;
// (*testout) << "p1, p2 = " << p1 << p2 << ", newp = " << p2d << endl;
newp = Point3d (p2d(0), p2d(1), 0);
newgi.edgenr = ap1.edgenr;
newgi.dist = newdist;
};
Vec<3> Refinement2d :: GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
Vec<2> t2d = geometry.GetSplines().Get(ap1.edgenr) -> GetTangent(ap1.dist);
return Vec<3> (t2d(0), t2d(1), 0);
}
Vec<3> Refinement2d :: GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const
{
return Vec<3> (0,0,1);
}
void Refinement2d :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & /* gi */) const
{
p(2) = 0;
}
void Refinement2d :: ProjectToEdge (Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & egi) const
{
Point<2> p2d (p(0), p(1)), pp;
double t;
geometry.GetSplines().Get(egi.edgenr) -> Project (p2d, pp, t);
p = Point<3> (pp(0), pp(1), 0);
}
}

View File

@ -1,52 +0,0 @@
#ifndef FILE_GEOM2DMESH
#define FILE_GEOM2DMESH
/**************************************************************************/
/* File: geom2dmesh.hh */
/* Author: Joachim Schoeberl */
/* Date: 22. Jan. 01 */
/**************************************************************************/
namespace netgen
{
class Refinement2d : public Refinement
{
const class SplineGeometry2d & geometry;
public:
Refinement2d (const class SplineGeometry2d & ageometry);
virtual ~Refinement2d ();
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
virtual Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const override;
virtual Vec<3> GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & /* gi */) const override;
virtual void ProjectToEdge (Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & egi) const override;
};
}
#endif

View File

@ -20,6 +20,76 @@ namespace netgen
delete [] materials[i];
}
void SplineGeometry2d :: PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
Point<2> p2d;
double newdist;
auto spline = GetSplines().Get(ap1.edgenr);
if( (ap1.dist == 0.0) && (ap2.dist == 0.0) )
{
// used for manually generated meshes
const SplineSeg3<2> * ss3;
const LineSeg<2> * ls;
auto ext = dynamic_cast<const SplineSegExt *>(spline);
if(ext)
{
ss3 = dynamic_cast<const SplineSeg3<2> *>(ext->seg);
ls = dynamic_cast<const LineSeg<2> *>(ext->seg);
}
else
{
ss3 = dynamic_cast<const SplineSeg3<2> *>(spline);
ls = dynamic_cast<const LineSeg<2> *>(spline);
}
Point<2> p12d(p1(0),p1(1)), p22d(p2(0),p2(1));
Point<2> p1_proj(0.0,0.0), p2_proj(0.0,0.0);
double t1_proj = 0.0;
double t2_proj = 0.0;
if(ss3)
{
ss3->Project(p12d,p1_proj,t1_proj);
ss3->Project(p22d,p2_proj,t2_proj);
}
else if(ls)
{
ls->Project(p12d,p1_proj,t1_proj);
ls->Project(p22d,p2_proj,t2_proj);
}
p2d = spline->GetPoint (((1-secpoint)*t1_proj+secpoint*t2_proj));
newdist = (1-secpoint)*t1_proj+secpoint*t2_proj;
}
else
{
p2d = spline->GetPoint (((1-secpoint)*ap1.dist+secpoint*ap2.dist));
newdist = (1-secpoint)*ap1.dist+secpoint*ap2.dist;
}
// (*testout) << "refine 2d line, ap1.dist, ap2.dist = " << ap1.dist << ", " << ap2.dist << endl;
// (*testout) << "p1, p2 = " << p1 << p2 << ", newp = " << p2d << endl;
newp = Point3d (p2d(0), p2d(1), 0);
newgi.edgenr = ap1.edgenr;
newgi.dist = newdist;
};
Vec<3> SplineGeometry2d :: GetTangent(const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
Vec<2> t2d = GetSplines().Get(ap1.edgenr) -> GetTangent(ap1.dist);
return Vec<3> (t2d(0), t2d(1), 0);
}
Vec<3> SplineGeometry2d :: GetNormal(int surfi1, const Point<3> & p,
const PointGeomInfo* gi) const
{
return Vec<3> (0,0,1);
}
void SplineGeometry2d :: Load (const char * filename)
{
@ -992,14 +1062,6 @@ namespace netgen
return 0;
}
Refinement & SplineGeometry2d :: GetRefinement () const
{
return * new Refinement2d (*this);
}
class SplineGeometryRegister : public GeometryRegister
{
public:

View File

@ -13,7 +13,6 @@
// #include "../gprim/spline.hpp"
// #include "../gprim/splinegeometry.hpp"
#include "geom2dmesh.hpp"
namespace netgen
{
@ -151,12 +150,40 @@ namespace netgen
void TestComment ( ifstream & infile ) ;
void DoArchive(Archive& ar)
void DoArchive(Archive& ar) override
{
SplineGeometry<2>::DoArchive(ar);
ar & materials & maxh & quadmeshing & tensormeshing & layer & bcnames & elto0;
}
bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override
{
p(2) = 0.0;
return true;
}
void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override
{
newp = p1+secpoint*(p2-p1);
newgi.trignum = 1;
}
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const override;
Vec<3> GetNormal(int surfi1, const Point<3> & p,
const PointGeomInfo* gi) const override;
const SplineSegExt & GetSpline (const int i) const
{
return dynamic_cast<const SplineSegExt&> (*splines[i]);
@ -168,7 +195,7 @@ namespace netgen
}
DLL_HEADER virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
DLL_HEADER int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam) override;
void PartitionBoundary (MeshingParameters & mp, double h, Mesh & mesh2d);
@ -205,9 +232,6 @@ namespace netgen
int AddBCName (string name);
string * BCNamePtr ( const int bcnr );
DLL_HEADER virtual Refinement & GetRefinement () const;
};
}

View File

@ -904,6 +904,14 @@ public:
GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
}
void GetIntersecting(const Point<dim> & pmin,
const Point<dim> & pmax,
Array<T> & pis) const
{
pis.SetSize0();
GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
}
void Insert (const Box<dim> & box, T pi)
{
Insert (box.PMin(), box.PMax(), pi);

View File

@ -19,6 +19,7 @@
#include <thread>
#include <mutex>
#include <atomic>
#include <optional>
#include <new>
#include <string>

View File

@ -37,5 +37,6 @@ install(FILES
localh.hpp meshclass.hpp meshfunc.hpp meshing2.hpp meshing3.hpp
meshing.hpp meshtool.hpp meshtype.hpp msghandler.hpp paralleltop.hpp
ruler2.hpp ruler3.hpp specials.hpp topology.hpp validate.hpp
python_mesh.hpp
DESTINATION ${NG_INSTALL_DIR_INCLUDE}/meshing COMPONENT netgen_devel
)

View File

@ -281,7 +281,7 @@ namespace netgen
NgArray<INDEX> & lindex,
double xh)
{
// static Timer timer("adfront2::GetLocals"); RegionTimer reg (timer);
static Timer timer("adfront2::GetLocals"); RegionTimer reg (timer);
int pstind;
Point<3> midp, p0;
@ -292,7 +292,7 @@ namespace netgen
loclines.Append(lines[baselineindex].L());
lindex.Append(baselineindex);
NgArrayMem<int, 1000> nearlines(0);
ArrayMem<int, 1000> nearlines(0);
NgArrayMem<int, 1000> nearpoints(0);
// dominating costs !!
@ -300,13 +300,14 @@ namespace netgen
p0 + Vec3d(xh, xh, xh),
nearlines);
pointsearchtree.GetIntersecting (p0 - Vec3d(xh, xh, xh),
// only special points that are not in adfront,
// other points are from linesearchtree
cpointsearchtree.GetIntersecting(p0 - Vec3d(xh, xh, xh),
p0 + Vec3d(xh, xh, xh),
nearpoints);
for (int ii = 0; ii < nearlines.Size(); ii++)
for(auto i : nearlines)
{
int i = nearlines[ii];
if (lines[i].Valid() && i != baselineindex)
{
loclines.Append(lines[i].L());
@ -317,38 +318,37 @@ namespace netgen
// static NgArray<int> invpindex;
invpindex.SetSize (points.Size());
// invpindex = -1;
for (int i = 0; i < nearpoints.Size(); i++)
invpindex[nearpoints[i]] = -1;
for(auto pi : nearpoints)
invpindex[pi] = -1;
for (int i = 0; i < loclines.Size(); i++)
for(const auto& li : loclines)
{
invpindex[loclines[i].I1()] = 0;
invpindex[loclines[i].I2()] = 0;
invpindex[li.I1()] = 0;
invpindex[li.I2()] = 0;
}
for (int i = 0; i < loclines.Size(); i++)
for(auto& line : loclines)
{
for (int j = 0; j < 2; j++)
for(auto i : Range(2))
{
int pi = loclines[i][j];
auto& pi = line[i];
if (invpindex[pi] == 0)
{
pindex.Append (pi);
invpindex[pi] = pindex.Size();
locpoints.Append (points[pi].P());
loclines[i][j] = locpoints.Size();
pi = locpoints.Size();
}
else
loclines[i][j] = invpindex[pi];
pi = invpindex[pi];
}
}
// double xh2 = xh*xh;
for (int ii = 0; ii < nearpoints.Size(); ii++)
for(auto i : nearpoints)
{
int i = nearpoints[ii];
if (points[i].Valid() &&
points[i].OnSurface() &&
// Dist2 (points.Get(i).P(), p0) <= xh2 &&

View File

@ -10,6 +10,499 @@ namespace netgen
GeometryRegister :: ~GeometryRegister()
{ ; }
void GeometryFace :: RestrictHTrig(Mesh& mesh,
const PointGeomInfo& gi0,
const PointGeomInfo& gi1,
const PointGeomInfo& gi2,
const MeshingParameters& mparam,
int depth, double h) const
{
auto p0 = GetPoint(gi0);
auto p1 = GetPoint(gi1);
auto p2 = GetPoint(gi2);
auto longest = (p0-p1).Length();
int cutedge = 2;
if(auto len = (p0-p2).Length(); len > longest)
{
longest = len;
cutedge = 1;
}
if(auto len = (p1-p2).Length(); len > longest)
{
longest = len;
cutedge = 0;
}
PointGeomInfo gi_mid;
gi_mid.u = (gi0.u + gi1.u + gi2.u)/3;
gi_mid.v = (gi0.v + gi1.v + gi2.v)/3;
if(depth % 3 == 0)
{
double curvature = 0.;
curvature = max({curvature, GetCurvature(gi_mid),
GetCurvature(gi0), GetCurvature(gi1),
GetCurvature(gi2)});
if(curvature < 1e-3)
return;
double kappa = curvature * mparam.curvaturesafety;
h = mparam.maxh * kappa < 1 ? mparam.maxh : 1./kappa;
if(h < 1e-4 * longest)
return;
}
if(h < longest && depth < 10)
{
if(cutedge == 0)
{
PointGeomInfo gi_m;
gi_m.u = 0.5 * (gi1.u + gi2.u);
gi_m.v = 0.5 * (gi1.v + gi2.v);
RestrictHTrig(mesh, gi_m, gi2, gi0, mparam, depth+1, h);
RestrictHTrig(mesh, gi_m, gi0, gi1, mparam, depth+1, h);
}
else if(cutedge == 1)
{
PointGeomInfo gi_m;
gi_m.u = 0.5 * (gi0.u + gi2.u);
gi_m.v = 0.5 * (gi0.v + gi2.v);
RestrictHTrig(mesh, gi_m, gi1, gi2, mparam, depth+1, h);
RestrictHTrig(mesh, gi_m, gi0, gi1, mparam, depth+1, h);
}
else if(cutedge == 2)
{
PointGeomInfo gi_m;
gi_m.u = 0.5 * (gi0.u + gi1.u);
gi_m.v = 0.5 * (gi0.v + gi1.v);
RestrictHTrig(mesh, gi_m, gi1, gi2, mparam, depth+1, h);
RestrictHTrig(mesh, gi_m, gi2, gi0, mparam, depth+1, h);
}
}
else
{
auto pmid = GetPoint(gi_mid);
for(const auto& p : {p0, p1, p2, pmid})
mesh.RestrictLocalH(p, h);
}
}
struct Line
{
Point<3> p0, p1;
inline double Length() const { return (p1-p0).Length(); }
inline double Dist(const Line& other) const
{
Vec<3> n = p1-p0;
Vec<3> q = other.p1-other.p0;
double nq = n*q;
Point<3> p = p0 + 0.5*n;
double lambda = (p-other.p0)*n / (nq + 1e-10);
if (lambda >= 0 && lambda <= 1)
return (p-other.p0-lambda*q).Length();
return 1e99;
}
};
void NetgenGeometry :: Analyse(Mesh& mesh,
const MeshingParameters& mparam) const
{
static Timer t1("SetLocalMeshsize"); RegionTimer regt(t1);
mesh.SetGlobalH(mparam.maxh);
mesh.SetMinimalH(mparam.minh);
mesh.SetLocalH(bounding_box.PMin(), bounding_box.PMax(),
mparam.grading);
// only set meshsize for edges longer than this
double mincurvelength = 1e-3 * bounding_box.Diam();
if(mparam.uselocalh)
{
double eps = 1e-10 * bounding_box.Diam();
const char* savetask = multithread.task;
multithread.task = "Analyse Edges";
// restrict meshsize on edges
for(auto i : Range(edges))
{
multithread.percent = 100. * i/edges.Size();
const auto & edge = edges[i];
auto length = edge->GetLength();
// skip very short edges
if(length < mincurvelength)
continue;
static constexpr int npts = 20;
// restrict mesh size based on edge length
for(auto i : Range(npts+1))
mesh.RestrictLocalH(edge->GetPoint(double(i)/npts), length/mparam.segmentsperedge);
// restrict mesh size based on edge curvature
double t = 0.;
auto p_old = edge->GetPoint(t);
while(t < 1.-eps)
{
t += edge->CalcStep(t, 1./mparam.curvaturesafety);
if(t < 1.)
{
auto p = edge->GetPoint(t);
auto dist = (p-p_old).Length();
mesh.RestrictLocalH(p, dist);
p_old = p;
}
}
}
multithread.task = "Analyse Faces";
// restrict meshsize on faces
for(auto i : Range(faces))
{
multithread.percent = 100. * i/faces.Size();
const auto& face = faces[i];
face->RestrictH(mesh, mparam);
}
if(mparam.closeedgefac.has_value())
{
multithread.task = "Analyse close edges";
constexpr int sections = 100;
Array<Line> lines;
lines.SetAllocSize(sections*edges.Size());
BoxTree<3> searchtree(bounding_box.PMin(),
bounding_box.PMax());
for(const auto& edge : edges)
{
if(edge->GetLength() < eps)
continue;
double t = 0.;
auto p_old = edge->GetPoint(t);
auto t_old = edge->GetTangent(t);
t_old.Normalize();
for(auto i : IntRange(1, sections+1))
{
t = double(i)/sections;
auto p_new = edge->GetPoint(t);
auto t_new = edge->GetTangent(t);
t_new.Normalize();
auto cosalpha = fabs(t_old * t_new);
if((i == sections) || (cosalpha < cos(10./180 * M_PI)))
{
auto index = lines.Append({p_old, p_new});
searchtree.Insert(p_old, p_new, index);
p_old = p_new;
t_old = t_new;
}
}
}
Array<int> linenums;
for(auto i : Range(lines))
{
const auto& line = lines[i];
if(line.Length() < eps) continue;
multithread.percent = 100.*i/lines.Size();
Box<3> box;
box.Set(line.p0);
box.Add(line.p1);
// box.Increase(max2(mesh.GetH(line.p0), mesh.GetH(line.p1)));
box.Increase(line.Length());
double mindist = 1e99;
linenums.SetSize0();
searchtree.GetIntersecting(box.PMin(), box.PMax(),
linenums);
for(auto num : linenums)
{
if(i == num) continue;
const auto & other = lines[num];
if((line.p0 - other.p0).Length2() < eps ||
(line.p0 - other.p1).Length2() < eps ||
(line.p1 - other.p0).Length2() < eps ||
(line.p1 - other.p1).Length2() < eps)
continue;
mindist = min2(mindist, line.Dist(other));
}
if(mindist == 1e99) continue;
mindist /= *mparam.closeedgefac + 1e-10;
if(mindist < 1e-3 * bounding_box.Diam())
{
(*testout) << "extremely small local h: " << mindist
<< " --> setting to " << 1e-3 * bounding_box.Diam() << endl;
(*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl
;
mindist = 1e-3 * bounding_box.Diam();
}
mesh.RestrictLocalHLine(line.p0, line.p1, mindist);
}
}
multithread.task = savetask;
}
for(const auto& mspnt : mparam.meshsize_points)
mesh.RestrictLocalH(mspnt.pnt, mspnt.h);
mesh.LoadLocalMeshSize(mparam.meshsizefilename);
}
void NetgenGeometry :: FindEdges(Mesh& mesh,
const MeshingParameters& mparam) const
{
static Timer t1("MeshEdges"); RegionTimer regt(t1);
static Timer tdivide("Divide Edges");
static Timer tdivedgesections("Divide edge sections");
const char* savetask = multithread.task;
multithread.task = "Mesh Edges";
// create face descriptors and set bc names
mesh.SetNBCNames(faces.Size());
for(auto i : Range(faces.Size()))
{
mesh.SetBCName(i, faces[i]->GetName());
// todo find attached solids
FaceDescriptor fd(i+1, 1, 0, i+1);
fd.SetBCName(mesh.GetBCNamePtr(i));
mesh.AddFaceDescriptor(fd);
}
std::map<size_t, PointIndex> vert2meshpt;
for(auto i : Range(vertices))
{
const auto& vert = *vertices[i];
MeshPoint mp(vert.GetPoint());
vert2meshpt[vert.GetHash()] = mesh.AddPoint(mp);
}
size_t segnr = 0;
for(auto facenr : Range(faces.Size()))
{
const auto& face = *faces[facenr];
for(auto facebndnr : Range(face.GetNBoundaries()))
{
auto boundary = face.GetBoundary(facebndnr);
for(auto enr : Range(boundary))
{
multithread.percent = 100. * ((double(enr)/boundary.Size() + facebndnr)/face.GetNBoundaries() + facenr)/faces.Size();
const auto& oriented_edge = *boundary[enr];
auto edgenr = GetEdgeIndex(oriented_edge);
const auto& edge = edges[edgenr];
PointIndex startp, endp;
// throws if points are not found
startp = vert2meshpt.at(edge->GetStartVertex().GetHash());
endp = vert2meshpt.at(edge->GetEndVertex().GetHash());
// ignore collapsed edges
if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam())
continue;
Array<MeshPoint> mps;
Array<double> params;
// -------------------- DivideEdge -----------------
static constexpr size_t divide_edge_sections = 1000;
tdivide.Start();
double hvalue[divide_edge_sections+1];
hvalue[0] = 0;
Point<3> old_pt = edge->GetPoint(0.);
// calc local h for edge
tdivedgesections.Start();
for(auto i : Range(divide_edge_sections))
{
auto pt = edge->GetPoint(double(i+1)/divide_edge_sections);
hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt) * (pt-old_pt).Length();
old_pt = pt;
}
int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
tdivedgesections.Stop();
mps.SetSize(nsubedges-1);
params.SetSize(nsubedges+1);
int i = 1;
int i1 = 0;
do
{
if (hvalue[i1]/hvalue[divide_edge_sections]*nsubedges >= i)
{
params[i] = (double(i1)/divide_edge_sections);
mps[i-1] = MeshPoint(edge->GetPoint(params[i]));
i++;
}
i1++;
if (i1 > divide_edge_sections)
{
nsubedges = i;
mps.SetSize(nsubedges-1);
params.SetSize(nsubedges+1);
cout << "divide edge: local h too small" << endl;
}
} while(i < nsubedges);
params[0] = 0.;
params[nsubedges] = 1.;
if(params[nsubedges] <= params[nsubedges-1])
{
cout << "CORRECTED" << endl;
mps.SetSize (nsubedges-2);
params.SetSize (nsubedges);
params[nsubedges-1] = 1.;
}
tdivide.Stop();
// ----------- Add Points to mesh and create segments -----
Array<PointIndex> pnums(mps.Size() + 2);
pnums[0] = startp;
pnums[mps.Size()+1] = endp;
double eps = bounding_box.Diam() * 1e-8;
for(auto i : Range(mps))
{
bool exists = false;
for(auto pi : Range(mesh.Points()))
{
if((mesh[pi] - mps[i]).Length() < eps)
{
exists = true;
pnums[i+1] = pi;
break;
}
}
if(!exists)
pnums[i+1] = mesh.AddPoint(mps[i]);
}
for(auto i : Range(pnums.Size()-1))
{
segnr++;
Segment seg;
seg[0] = pnums[i];
seg[1] = pnums[i+1];
seg.edgenr = segnr;
seg.epgeominfo[0].dist = params[i];
seg.epgeominfo[1].dist = params[i+1];
seg.epgeominfo[0].edgenr = edgenr;
seg.epgeominfo[1].edgenr = edgenr;
seg.si = facenr+1;
seg.surfnr1 = facenr+1;
// TODO: implement functionality to transfer edge parameter t to face parameters u,v
for(auto j : Range(2))
face.CalcEdgePointGI(*edge, params[i+j],
seg.epgeominfo[j]);
if(!oriented_edge.OrientedLikeGlobal())
{
swap (seg[0], seg[1]);
swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
}
mesh.AddSegment(seg);
}
}
}
}
mesh.CalcSurfacesOfNode();
multithread.task = savetask;
}
void NetgenGeometry :: MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) const
{
static Timer t1("Surface Meshing"); RegionTimer regt(t1);
const char* savetask = multithread.task;
multithread.task = "Mesh Surface";
Array<int, PointIndex> glob2loc(mesh.GetNP());
for(auto k : Range(faces))
{
multithread.percent = 100. * k/faces.Size();
const auto& face = *faces[k];
auto bb = face.GetBoundingBox();
bb.Increase(bb.Diam()/10);
Meshing2 meshing(*this, mparam, bb);
glob2loc = 0;
int cntp = 0;
for(auto& seg : mesh.LineSegments())
{
if(seg.si == k+1)
{
for(auto j : Range(2))
{
auto pi = seg[j];
if(glob2loc[pi] == 0)
{
meshing.AddPoint(mesh[pi], pi);
cntp++;
glob2loc[pi] = cntp;
}
}
}
}
for(auto & seg : mesh.LineSegments())
{
if(seg.si == k+1)
{
PointGeomInfo gi0, gi1;
gi0.trignum = gi1.trignum = k+1;
gi0.u = seg.epgeominfo[0].u;
gi0.v = seg.epgeominfo[0].v;
gi1.u = seg.epgeominfo[1].u;
gi1.v = seg.epgeominfo[1].v;
meshing.AddBoundaryElement(glob2loc[seg[0]],
glob2loc[seg[1]],
gi0, gi1);
}
}
// TODO Set max area 2* area of face
auto noldsurfels = mesh.GetNSE();
static Timer t("GenerateMesh"); RegionTimer reg(t);
MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1);
for(auto i : Range(noldsurfels, mesh.GetNSE()))
{
mesh.SurfaceElements()[i].SetIndex(k+1);
}
}
multithread.task = savetask;
}
void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const
{
const auto savetask = multithread.task;
multithread.task = "Optimizing surface";
static Timer timer_opt2d("Optimization 2D");
RegionTimer reg(timer_opt2d);
auto meshopt = MeshOptimize2d(mesh);
for(auto i : Range(mparam.optsteps2d))
{
PrintMessage(3, "Optimization step ", i);
int innerstep = 0;
for(auto optstep : mparam.optimize2d)
{
multithread.percent = 100. * (double(innerstep++)/mparam.optimize2d.size() + i)/mparam.optsteps2d;
switch(optstep)
{
case 's':
meshopt.EdgeSwapping(0);
break;
case 'S':
meshopt.EdgeSwapping(1);
break;
case 'm':
meshopt.ImproveMesh(mparam);
break;
case 'c':
meshopt.CombineImprove();
break;
}
}
}
mesh.CalcSurfacesOfNode();
mesh.Compress();
multithread.task = savetask;
}
shared_ptr<NetgenGeometry> GeometryRegisterArray :: LoadFromMeshFile (istream & ist) const
{
@ -27,17 +520,48 @@ namespace netgen
int NetgenGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
{
if (!mesh) return 1;
multithread.percent = 0;
if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
if(mparam.perfstepsstart <= MESHCONST_ANALYSE)
{
if(!mesh)
mesh = make_shared<Mesh>();
mesh->geomtype = GetGeomType();
Analyse(*mesh, mparam);
}
if(multithread.terminate || mparam.perfstepsend <= MESHCONST_ANALYSE)
return 0;
if(mparam.perfstepsstart <= MESHCONST_MESHEDGES)
FindEdges(*mesh, mparam);
if(multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES)
return 0;
if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
{
MeshSurface(*mesh, mparam);
mesh->CalcSurfacesOfNode();
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHSURFACE)
return 0;
if (mparam.perfstepsstart <= MESHCONST_OPTSURFACE)
OptimizeSurface(*mesh, mparam);
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE)
return 0;
if(mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
{
multithread.task = "Volume meshing";
MESHING3_RESULT res =
MeshVolume (mparam, *mesh);
MESHING3_RESULT res = MeshVolume (mparam, *mesh);
if (res != MESHING3_OK) return 1;
if (multithread.terminate) return 0;
RemoveIllegalElements (*mesh);
@ -46,7 +570,6 @@ namespace netgen
MeshQuality3d (*mesh);
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHVOLUME)
return 0;
@ -58,17 +581,10 @@ namespace netgen
OptimizeVolume (mparam, *mesh);
if (multithread.terminate) return 0;
}
FinalizeMesh(*mesh);
return 0;
}
const Refinement & NetgenGeometry :: GetRefinement () const
{
return *new Refinement;;
}
void NetgenGeometry :: Save (string filename) const
{
throw NgException("Cannot save geometry - no geometry available");

View File

@ -12,19 +12,181 @@ struct Tcl_Interp;
namespace netgen
{
class GeometryVertex
{
public:
virtual ~GeometryVertex() {}
virtual Point<3> GetPoint() const = 0;
virtual size_t GetHash() const = 0;
};
class GeometryEdge
{
public:
virtual ~GeometryEdge() {}
virtual const GeometryVertex& GetStartVertex() const = 0;
virtual const GeometryVertex& GetEndVertex() const = 0;
virtual double GetLength() const = 0;
virtual Point<3> GetPoint(double t) const = 0;
// Calculate parameter step respecting edges sag value
virtual double CalcStep(double t, double sag) const = 0;
virtual bool OrientedLikeGlobal() const = 0;
virtual size_t GetHash() const = 0;
virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const = 0;
virtual void PointBetween(const Point<3>& p1,
const Point<3>& p2,
double secpoint,
const EdgePointGeomInfo& gi1,
const EdgePointGeomInfo& gi2,
Point<3>& newp,
EdgePointGeomInfo& newgi) const = 0;
virtual Vec<3> GetTangent(double t) const = 0;
};
class GeometryFace
{
public:
virtual ~GeometryFace() {}
virtual size_t GetNBoundaries() const = 0;
virtual Array<unique_ptr<GeometryEdge>> GetBoundary(size_t index) const = 0;
virtual string GetName() const { return "default"; }
virtual PointGeomInfo Project(Point<3>& p) const = 0;
// Project point using geo info. Fast if point is close to
// parametrization in geo info.
virtual bool ProjectPointGI(Point<3>& p, PointGeomInfo& gi) const =0;
virtual bool CalcPointGeomInfo(const Point<3>& p, PointGeomInfo& gi) const = 0;
virtual Point<3> GetPoint(const PointGeomInfo& gi) const = 0;
virtual void CalcEdgePointGI(const GeometryEdge& edge,
double t,
EdgePointGeomInfo& egi) const = 0;
virtual Box<3> GetBoundingBox() const = 0;
// Get curvature in point from local coordinates in PointGeomInfo
virtual double GetCurvature(const PointGeomInfo& gi) const = 0;
virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const = 0;
virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo* gi = nullptr) const = 0;
virtual void PointBetween(const Point<3>& p1,
const Point<3>& p2,
double secpoint,
const PointGeomInfo& gi1,
const PointGeomInfo& gi2,
Point<3>& newp,
PointGeomInfo& newgi) const = 0;
protected:
void RestrictHTrig(Mesh& mesh,
const PointGeomInfo& gi0,
const PointGeomInfo& gi1,
const PointGeomInfo& gi2,
const MeshingParameters& mparam,
int depth = 0, double h = 0.) const;
};
class DLL_HEADER NetgenGeometry
{
unique_ptr<Refinement> ref;
protected:
Array<unique_ptr<GeometryVertex>> vertices;
Array<unique_ptr<GeometryEdge>> edges;
Array<unique_ptr<GeometryFace>> faces;
Box<3> bounding_box;
public:
NetgenGeometry()
{
ref = make_unique<Refinement>(*this);
}
virtual ~NetgenGeometry () { ; }
virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
virtual const Refinement & GetRefinement () const;
virtual const Refinement & GetRefinement () const
{
return *ref;
}
virtual void DoArchive(Archive&)
{ throw NgException("DoArchive not implemented for " + Demangle(typeid(*this).name())); }
virtual Mesh::GEOM_TYPE GetGeomType() const { return Mesh::NO_GEOM; }
virtual void Analyse(Mesh& mesh,
const MeshingParameters& mparam) const;
virtual void FindEdges(Mesh& mesh, const MeshingParameters& mparam) const;
virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) const;
virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const;
virtual void FinalizeMesh(Mesh& mesh) const {}
virtual PointGeomInfo ProjectPoint (int surfind, Point<3> & p) const
{
return faces[surfind-1]->Project(p);
}
virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo* gi = nullptr) const
{
edges[gi->edgenr]->ProjectPoint(p, gi);
}
virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const
{
return faces[surfind-1]->CalcPointGeomInfo(p3, gi);
}
virtual bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
{
return faces[surfind-1]->ProjectPointGI(p, gi);
}
virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const
{ return faces[surfind-1]->GetNormal(p, gi); }
virtual void PointBetween (const Point<3> & p1,
const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp,
PointGeomInfo & newgi) const
{
if(faces.Size())
{
faces[surfi-1]->PointBetween(p1, p2, secpoint, gi1, gi2, newp, newgi);
return;
}
newp = p1 + secpoint * (p2-p1);
}
virtual void PointBetweenEdge(const Point<3> & p1,
const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp,
EdgePointGeomInfo & newgi) const
{
if(edges.Size())
{
edges[ap1.edgenr]->PointBetween(p1, p2, secpoint,
ap1, ap2, newp, newgi);
return;
}
newp = p1+secpoint*(p2-p1);
}
virtual Vec<3> GetTangent(const Point<3> & p, int surfi1,
int surfi2,
const EdgePointGeomInfo & egi) const
{
throw Exception("Base geometry get tangent called");
}
virtual size_t GetEdgeIndex(const GeometryEdge& edge) const
{
for(auto i : Range(edges))
if(edge.GetHash() == edges[i]->GetHash())
return i;
throw Exception("Couldn't find edge index");
}
virtual void Save (string filename) const;
virtual void SaveToMeshFile (ostream & /* ost */) const { ; }
};

View File

@ -1958,6 +1958,8 @@ namespace netgen
const NgArray< NgArray<int,PointIndex::BASE>* > & idmaps,
const string & refinfofile)
{
if (mesh.GetDimension() < 2)
throw Exception ("Mesh bisection is available in 2D and 3D");
// mtets.SetName ("bisection, tets");
// mprisms.SetName ("bisection, prisms");
// mtris.SetName ("bisection, trigs");
@ -3428,7 +3430,7 @@ namespace netgen
PointGeomInfo npgi;
if (mesh[newp].Type() != EDGEPOINT)
PointBetween (mesh.Point (oldpi1), mesh.Point (oldpi2),
geo.PointBetween (mesh.Point (oldpi1), mesh.Point (oldpi2),
0.5, si,
oldtri.pgeominfo[(oldtri.markededge+1)%3],
oldtri.pgeominfo[(oldtri.markededge+2)%3],
@ -3506,28 +3508,16 @@ namespace netgen
PointGeomInfo npgi1, npgi2;
int si = mesh.GetFaceDescriptor (oldquad.surfid).SurfNr();
// geom->GetSurface(si)->Project (mesh.Point(newp1));
// geom->GetSurface(si)->Project (mesh.Point(newp2));
// (*testout)
// cerr << "project point 1 " << newp1 << " old: " << mesh.Point(newp1);
PointBetween (mesh.Point (edge1.I1()), mesh.Point (edge1.I2()),
geo.PointBetween(mesh.Point (edge1.I1()), mesh.Point (edge1.I2()),
0.5, si,
pgi11,
pgi12,
mesh.Point (newp1), npgi1);
// (*testout)
// cerr << " new: " << mesh.Point(newp1) << endl;
// cerr << "project point 2 " << newp2 << " old: " << mesh.Point(newp2);
PointBetween (mesh.Point (edge2.I1()), mesh.Point (edge2.I2()),
geo.PointBetween (mesh.Point (edge2.I1()), mesh.Point (edge2.I2()),
0.5, si,
pgi21,
pgi22,
mesh.Point (newp2), npgi2);
// cerr << " new: " << mesh.Point(newp2) << endl;
BTBisectQuad (oldquad, newp1, npgi1, newp2, npgi2,
newquad1, newquad2);
@ -3563,16 +3553,10 @@ namespace netgen
EdgePointGeomInfo newepgi;
//
// cerr << "move edgepoint " << newpi << " from " << mesh.Point(newpi);
PointBetween (mesh.Point (seg[0]), mesh.Point (seg[1]),
geo.PointBetweenEdge(mesh.Point (seg[0]), mesh.Point (seg[1]),
0.5, seg.surfnr1, seg.surfnr2,
seg.epgeominfo[0], seg.epgeominfo[1],
mesh.Point (newpi), newepgi);
// cerr << " to " << mesh.Point (newpi) << endl;
nseg1.epgeominfo[1] = newepgi;
nseg2.epgeominfo[0] = newepgi;
@ -4139,62 +4123,4 @@ namespace netgen
refine_hp = 0;
refine_p = 0;
}
Refinement :: Refinement ()
{
optimizer2d = NULL;
}
Refinement :: ~Refinement ()
{
;
}
void Refinement :: PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
}
void Refinement :: PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
//cout << "base class edge point between" << endl;
newp = p1+secpoint*(p2-p1);
}
Vec<3> Refinement :: GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & ap1) const
{
cerr << "Refinement::GetTangent not overloaded" << endl;
return Vec<3> (0,0,0);
}
Vec<3> Refinement :: GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const
{
cerr << "Refinement::GetNormal not overloaded" << endl;
return Vec<3> (0,0,0);
}
void Refinement :: ProjectToSurface (Point<3> & p, int surfi) const
{
if (printmessage_importance>0)
cerr << "Refinement :: ProjectToSurface ERROR: no geometry set" << endl;
};
void Refinement :: ProjectToEdge (Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const
{
cerr << "Refinement::ProjectToEdge not overloaded" << endl;
}
}

View File

@ -38,11 +38,11 @@ DLL_HEADER extern void ZRefinement (Mesh &, const class NetgenGeometry *,
class DLL_HEADER Refinement
{
MeshOptimize2d * optimizer2d;
const NetgenGeometry& geo;
public:
Refinement ();
virtual ~Refinement ();
Refinement (const NetgenGeometry& ageo) : geo(ageo) {}
virtual ~Refinement () {}
void Refine (Mesh & mesh) const;
void Refine (Mesh & mesh);
@ -51,49 +51,10 @@ public:
void MakeSecondOrder (Mesh & mesh) const;
void MakeSecondOrder (Mesh & mesh);
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const;
virtual Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2,
const EdgePointGeomInfo & egi) const;
virtual Vec<3> GetNormal (const Point<3> & p, int surfi1,
const PointGeomInfo & gi) const;
virtual void ProjectToSurface (Point<3> & p, int surfi) const;
virtual void ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & /* gi */) const
{
ProjectToSurface (p, surfi);
}
virtual void ProjectToEdge (Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const;
void ValidateSecondOrder (Mesh & mesh);
void ValidateRefinedMesh (Mesh & mesh,
NgArray<INDEX_2> & parents);
MeshOptimize2d * Get2dOptimizer(void) const
{
return optimizer2d;
}
void Set2dOptimizer(MeshOptimize2d * opti)
{
optimizer2d = opti;
}
virtual void LocalizeEdgePoints(Mesh & /* mesh */) const {;}
};

View File

@ -539,7 +539,7 @@ namespace netgen
CurvedElements :: CurvedElements (const Mesh & amesh)
: mesh (amesh)
: mesh(amesh)
{
order = 1;
rational = 0;
@ -555,6 +555,7 @@ namespace netgen
void CurvedElements :: BuildCurvedElements(const Refinement * ref, int aorder,
bool arational)
{
auto & geo = *mesh.GetGeometry();
ishighorder = 0;
order = 1;
@ -838,8 +839,8 @@ namespace netgen
{
Point<3> pm = Center (p1, p2);
Vec<3> n1 = ref -> GetNormal (p1, surfnr[e], gi0[e]);
Vec<3> n2 = ref -> GetNormal (p2, surfnr[e], gi1[e]);
Vec<3> n1 = geo.GetNormal (surfnr[e], p1, &gi0[e]);
Vec<3> n2 = geo.GetNormal (surfnr[e], p2, &gi1[e]);
// p3 = pm + alpha1 n1 + alpha2 n2
@ -876,7 +877,7 @@ namespace netgen
Vec<3> v05 = 0.25 * Vec<3> (p1) + 0.5*w* Vec<3>(p3) + 0.25 * Vec<3> (p2);
v05 /= 1 + (w-1) * 0.5;
Point<3> p05 (v05), pp05(v05);
ref -> ProjectToSurface (pp05, surfnr[e], gi0[e]);
geo.ProjectPointGI(surfnr[e], pp05, gi0[e]);
double d = Dist (pp05, p05);
if (d < dold)
@ -911,14 +912,14 @@ namespace netgen
if (swap)
{
p = p1 + xi[j] * (p2-p1);
ref -> PointBetween (p1, p2, xi[j],
geo.PointBetween (p1, p2, xi[j],
surfnr[e], gi0[e], gi1[e],
pp, ppgi);
}
else
{
p = p2 + xi[j] * (p1-p2);
ref -> PointBetween (p2, p1, xi[j],
geo.PointBetween (p2, p1, xi[j],
surfnr[e], gi1[e], gi0[e],
pp, ppgi);
}
@ -1053,9 +1054,9 @@ namespace netgen
if (rational)
{
Vec<3> tau1 = ref -> GetTangent (p1, edge_surfnr2[edgenr], edge_surfnr1[edgenr],
Vec<3> tau1 = geo.GetTangent(p1, edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi0[edgenr]);
Vec<3> tau2 = ref -> GetTangent (p2, edge_surfnr2[edgenr], edge_surfnr1[edgenr],
Vec<3> tau2 = geo.GetTangent(p2, edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi1[edgenr]);
// p1 + alpha1 tau1 = p2 + alpha2 tau2;
@ -1082,8 +1083,8 @@ namespace netgen
Vec<3> v05 = 0.25 * Vec<3> (p1) + 0.5*w* Vec<3>(p3) + 0.25 * Vec<3> (p2);
v05 /= 1 + (w-1) * 0.5;
Point<3> p05 (v05), pp05(v05);
ref -> ProjectToEdge (pp05, edge_surfnr1[edgenr], edge_surfnr2[edgenr],
edge_gi0[edgenr]);
geo.ProjectPointEdge(edge_surfnr1[edgenr], edge_surfnr2[edgenr], pp05,
&edge_gi0[edgenr]);
double d = Dist (pp05, p05);
if (d < dold)
@ -1127,7 +1128,7 @@ namespace netgen
if (swap)
{
p = p1 + xi[j] * (p2-p1);
ref -> PointBetween (p1, p2, xi[j],
geo.PointBetweenEdge(p1, p2, xi[j],
edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi0[edgenr], edge_gi1[edgenr],
pp, ppgi);
@ -1135,7 +1136,7 @@ namespace netgen
else
{
p = p2 + xi[j] * (p1-p2);
ref -> PointBetween (p2, p1, xi[j],
geo.PointBetweenEdge(p2, p1, xi[j],
edge_surfnr2[edgenr], edge_surfnr1[edgenr],
edge_gi1[edgenr], edge_gi0[edgenr],
pp, ppgi);
@ -1302,10 +1303,10 @@ namespace netgen
SurfaceElementIndex sei = top.GetFace2SurfaceElement (f+1)-1;
if (sei != SurfaceElementIndex(-1)) {
PointGeomInfo gi = mesh[sei].GeomInfoPi(1);
ref -> ProjectToSurface (pp, surfnr[facenr], gi);
geo.ProjectPointGI(surfnr[facenr], pp, gi);
}
else
{ ref -> ProjectToSurface (pp, surfnr[facenr]); }
{ geo.ProjectPoint(surfnr[facenr], pp); }
Vec<3> dist = pp-xa[jj];
CalcTrigShape (order1, lami[fnums[1]]-lami[fnums[0]],

View File

@ -1,9 +1,289 @@
#include <mystdlib.h>
#include "meshing.hpp"
namespace netgen
{
template<int dim, typename T=INDEX, typename TSCAL=double>
class DelaunayTree
{
public:
// Number of entries per leaf
static constexpr int N = 100;
struct Node;
struct Leaf
{
Point<2*dim, TSCAL> p[N];
T index[N];
int n_elements;
int nr;
Leaf() : n_elements(0)
{ }
void Add( Array<Leaf*> &leaves, Array<T> &leaf_index, const Point<2*dim> &ap, T aindex )
{
p[n_elements] = ap;
index[n_elements] = aindex;
n_elements++;
if(leaf_index.Size()<aindex+1)
leaf_index.SetSize(aindex+1);
leaf_index[aindex] = nr;
}
};
struct Node
{
union
{
Node *children[2];
Leaf *leaf;
};
double sep;
int level;
Node()
: children{nullptr,nullptr}
{ }
~Node()
{ }
Leaf *GetLeaf() const
{
return children[1] ? nullptr : leaf;
}
};
private:
Node root;
Array<Leaf*> leaves;
Array<T> leaf_index;
Point<dim> global_min, global_max;
double tol;
size_t n_leaves;
size_t n_nodes;
BlockAllocator ball_nodes;
BlockAllocator ball_leaves;
public:
DelaunayTree (const Point<dim> & pmin, const Point<dim> & pmax)
: global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf))
{
root.leaf = (Leaf*) ball_leaves.Alloc(); new (root.leaf) Leaf();
root.leaf->nr = 0;
leaves.Append(root.leaf);
root.level = 0;
tol = 1e-7 * Dist(pmax, pmin);
}
DelaunayTree (const Box<dim> & box)
: DelaunayTree(box.PMin(), box.PMax())
{ }
size_t GetNLeaves()
{
return n_leaves;
}
size_t GetNNodes()
{
return n_nodes;
}
template<typename TFunc>
void GetFirstIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
TFunc func=[](auto pi){return false;}) const
{
// static Timer timer("DelaunayTree::GetIntersecting"); RegionTimer rt(timer);
// static Timer timer1("DelaunayTree::GetIntersecting-LinearSearch");
ArrayMem<const Node*, 100> stack;
ArrayMem<int, 100> dir_stack;
Point<2*dim> tpmin, tpmax;
for (size_t i : IntRange(dim))
{
tpmin(i) = global_min(i);
tpmax(i) = pmax(i)+tol;
tpmin(i+dim) = pmin(i)-tol;
tpmax(i+dim) = global_max(i);
}
stack.SetSize(0);
stack.Append(&root);
dir_stack.SetSize(0);
dir_stack.Append(0);
while(stack.Size())
{
const Node *node = stack.Last();
stack.DeleteLast();
int dir = dir_stack.Last();
dir_stack.DeleteLast();
if(Leaf *leaf = node->GetLeaf())
{
// RegionTimer rt1(timer1);
for (auto i : IntRange(leaf->n_elements))
{
bool intersect = true;
const auto p = leaf->p[i];
for (int d = 0; d < dim; d++)
if (p[d] > tpmax[d])
intersect = false;
for (int d = dim; d < 2*dim; d++)
if (p[d] < tpmin[d])
intersect = false;
if(intersect)
if(func(leaf->index[i])) return;
}
}
else
{
int newdir = dir+1;
if(newdir==2*dim) newdir = 0;
if (tpmin[dir] <= node->sep)
{
stack.Append(node->children[0]);
dir_stack.Append(newdir);
}
if (tpmax[dir] >= node->sep)
{
stack.Append(node->children[1]);
dir_stack.Append(newdir);
}
}
}
}
void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
NgArray<T> & pis) const
{
pis.SetSize(0);
GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
}
void Insert (const Box<dim> & box, T pi)
{
Insert (box.PMin(), box.PMax(), pi);
}
void Insert (const Point<dim> & pmin, const Point<dim> & pmax, T pi)
{
// static Timer timer("DelaunayTree::Insert"); RegionTimer rt(timer);
int dir = 0;
Point<2*dim> p;
for (auto i : IntRange(dim))
{
p(i) = pmin[i];
p(i+dim) = pmax[i];
}
Node * node = &root;
Leaf * leaf = node->GetLeaf();
// search correct leaf to add point
while(!leaf)
{
node = p[dir] < node->sep ? node->children[0] : node->children[1];
dir++;
if(dir==2*dim) dir = 0;
leaf = node->GetLeaf();
}
// add point to leaf
if(leaf->n_elements < N)
leaf->Add(leaves, leaf_index, p,pi);
else // assume leaf->n_elements == N
{
// add two new nodes and one new leaf
int n_elements = leaf->n_elements;
ArrayMem<TSCAL, N> coords(n_elements);
ArrayMem<int, N> order(n_elements);
// separate points in two halves, first sort all coordinates in direction dir
for (auto i : IntRange(n_elements))
{
order[i] = i;
coords[i] = leaf->p[i][dir];
}
QuickSortI(coords, order);
int isplit = N/2;
Leaf *leaf1 = (Leaf*) ball_leaves.Alloc(); new (leaf1) Leaf();
Leaf *leaf2 = (Leaf*) ball_leaves.Alloc(); new (leaf2) Leaf();
leaf1->nr = leaf->nr;
leaf2->nr = leaves.Size();
leaves.Append(leaf2);
leaves[leaf1->nr] = leaf1;
for (auto i : order.Range(isplit))
leaf1->Add(leaves, leaf_index, leaf->p[i], leaf->index[i] );
for (auto i : order.Range(isplit, N))
leaf2->Add(leaves, leaf_index, leaf->p[i], leaf->index[i] );
Node *node1 = (Node*) ball_nodes.Alloc(); new (node1) Node();
node1->leaf = leaf1;
node1->level = node->level+1;
Node *node2 = (Node*) ball_nodes.Alloc(); new (node2) Node();
node2->leaf = leaf2;
node2->level = node->level+1;
node->children[0] = node1;
node->children[1] = node2;
node->sep = 0.5 * (leaf->p[order[isplit-1]][dir] + leaf->p[order[isplit]][dir]);
// add new point to one of the new leaves
if (p[dir] < node->sep)
leaf1->Add( leaves, leaf_index, p, pi );
else
leaf2->Add( leaves, leaf_index, p, pi );
ball_leaves.Free(leaf);
n_leaves++;
n_nodes+=2;
}
}
void DeleteElement (T pi)
{
// static Timer timer("DelaunayTree::DeleteElement"); RegionTimer rt(timer);
Leaf *leaf = leaves[leaf_index[pi]];
leaf_index[pi] = -1;
auto & n_elements = leaf->n_elements;
auto & index = leaf->index;
auto & p = leaf->p;
for (auto i : IntRange(n_elements))
{
if(index[i] == pi)
{
n_elements--;
if(i!=n_elements)
{
index[i] = index[n_elements];
p[i] = p[n_elements];
}
return;
}
}
}
};
// typedef BoxTree<3> DTREE;
typedef DelaunayTree<3> DTREE;
static const int deltetfaces[][3] =
{ { 1, 2, 3 },
@ -227,14 +507,14 @@ namespace netgen
void AddDelaunayPoint (PointIndex newpi, const Point3d & newp,
NgArray<DelaunayTet> & tempels,
Mesh & mesh,
BoxTree<3> & tettree,
DTREE & tettree,
MeshNB & meshnb,
NgArray<Point<3> > & centers, NgArray<double> & radi2,
NgArray<int> & connected, NgArray<int> & treesearch,
NgArray<int> & freelist, SphereList & list,
IndexSet & insphere, IndexSet & closesphere)
{
static Timer t("Meshing3::AddDelaunayPoint"); RegionTimer reg(t);
static Timer t("Meshing3::AddDelaunayPoint");// RegionTimer reg(t);
// static Timer tsearch("addpoint, search");
// static Timer tinsert("addpoint, insert");
@ -635,7 +915,7 @@ namespace netgen
pmin2 = pmin2 + 0.1 * (pmin2 - pmax2);
pmax2 = pmax2 + 0.1 * (pmax2 - pmin2);
BoxTree<3> tettree(pmin2, pmax2);
DTREE tettree(pmin2, pmax2);
tempels.Append (startel);
@ -798,6 +1078,7 @@ namespace netgen
tempmesh.AddVolumeElement (el);
}
tempels.DeleteAll();
MeshQuality3d (tempmesh);
@ -852,6 +1133,7 @@ namespace netgen
MeshQuality3d (tempmesh);
tempels.SetSize(tempmesh.GetNE());
tempels.SetSize(0);
for (auto & el : tempmesh.VolumeElements())
tempels.Append (el);

View File

@ -118,10 +118,10 @@ inline int FindInnerPoint (POINTArray & points,
{
// const Element2d & el = faces[i];
// (*testout) << "el[" << i << "] = " << el << endl;
for (int j = 1; j <= 3; j++)
for (int j : Range(3))
{
double hi = Dist (points[faces[i].PNumMod(j)],
points[faces[i].PNumMod(j+1)]);
double hi = Dist (points[faces[i][j%3]],
points[faces[i][(j+1)%3]]);
if (hi > hmax) hmax = hi;
}
}

View File

@ -77,8 +77,6 @@ namespace netgen
NgArray<int> tets_in_qualclass;
mutex tcl_todo_mutex;
int h_argc = 0;

View File

@ -27,8 +27,6 @@ namespace netgen
// extern DLL_HEADER MeshingParameters mparam;
DLL_HEADER extern NgArray<int> tets_in_qualclass;
DLL_HEADER extern mutex tcl_todo_mutex;
class DLL_HEADER multithreadt

File diff suppressed because it is too large Load Diff

View File

@ -1,32 +1,105 @@
#ifndef FILE_IMPROVE2
#define FILE_IMPROVE2
template<typename TINDEX>
void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges )
{
static Timer tbuild_edges("Build edges"); RegionTimer reg(tbuild_edges);
static constexpr int tetedges[6][2] =
{ { 0, 1 }, { 0, 2 }, { 0, 3 },
{ 1, 2 }, { 1, 3 }, { 2, 3 } };
int ntasks = 2*ngcore::TaskManager::GetMaxThreads();
Array<Array<std::tuple<PointIndex,PointIndex>>> task_edges(ntasks);
ParallelFor(IntRange(ntasks), [&] (int ti)
{
auto myrange = mesh.Points().Range().Split(ti, ntasks);
ArrayMem<std::tuple<PointIndex,PointIndex>, 100> local_edges;
for (auto pi : myrange)
{
local_edges.SetSize(0);
for(auto ei : elementsonnode[pi])
{
const auto & elem = mesh[ei];
if (elem.IsDeleted()) continue;
for (int j = 0; j < 6; j++)
{
PointIndex pi0 = elem[tetedges[j][0]];
PointIndex pi1 = elem[tetedges[j][1]];
if (pi1 < pi0) Swap(pi0, pi1);
if(pi0==pi)
local_edges.Append(std::make_tuple(pi0, pi1));
}
}
QuickSort(local_edges);
auto edge_prev = std::make_tuple<PointIndex, PointIndex>(-1,-1);
for(auto edge : local_edges)
if(edge != edge_prev)
{
task_edges[ti].Append(edge);
edge_prev = edge;
}
}
}, ntasks);
int num_edges = 0;
for (auto & edg : task_edges)
num_edges += edg.Size();
edges.SetAllocSize(num_edges);
for (auto & edg : task_edges)
edges.Append(edg);
}
class Neighbour
{
int nr[3];
int orient[3];
public:
Neighbour () { ; }
void SetNr (int side, int anr) { nr[side] = anr; }
int GetNr (int side) { return nr[side]; }
void SetOrientation (int side, int aorient) { orient[side] = aorient; }
int GetOrientation (int side) { return orient[side]; }
};
///
class MeshOptimize2d
{
int faceindex;
int improveedges;
double metricweight;
int writestatus;
int faceindex = 0;
int improveedges = 0;
double metricweight = 0.;
int writestatus = 1;
Mesh& mesh;
const NetgenGeometry& geo;
public:
///
MeshOptimize2d ();
MeshOptimize2d(Mesh& amesh) : mesh(amesh), geo(*mesh.GetGeometry())
{}
virtual ~MeshOptimize2d() { ; }
///
void ImproveMesh (Mesh & mesh2d, const MeshingParameters & mp);
void ImproveMeshJacobian (Mesh & mesh2d, const MeshingParameters & mp);
void ImproveVolumeMesh (Mesh & mesh);
void ImproveMesh (const MeshingParameters & mp);
void ImproveMeshJacobian (const MeshingParameters & mp);
void ImproveVolumeMesh ();
void ProjectBoundaryPoints(NgArray<int> & surfaceindex,
const NgArray<Point<3>* > & from, NgArray<Point<3>* > & dest);
void EdgeSwapping (Mesh & mesh, int usemetric);
void CombineImprove (Mesh & mesh);
void SplitImprove (Mesh & mesh);
bool EdgeSwapping (const int usemetric, Array<Neighbour> &neighbors, Array<bool> &swapped,
const SurfaceElementIndex t1, const int edge, const int t, Array<int,PointIndex> &pdef, const bool check_only=false);
void EdgeSwapping (int usemetric);
void CombineImprove ();
void SplitImprove ();
void GenericImprove (Mesh & mesh);
void GenericImprove ();
void SetFaceIndex (int fi) { faceindex = fi; }
@ -35,31 +108,9 @@ public:
void SetWriteStatus (int ws) { writestatus = ws; }
///
virtual void SelectSurfaceOfPoint (const Point<3> & p,
const PointGeomInfo & gi);
///
virtual void ProjectPoint (INDEX /* surfind */, Point<3> & /* p */) const { };
/// project point, use gi as initial value, and compute new gi
virtual int ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
{ ProjectPoint (surfind, p); return CalcPointGeomInfo (surfind, gi, p); }
///
virtual void ProjectPoint2 (INDEX /* surfind */, INDEX /* surfind2 */, Point<3> & /* p */) const { };
/// liefert zu einem 3d-Punkt die geominfo (Dreieck) und liefert 1, wenn erfolgreich,
/// 0, wenn nicht (Punkt ausserhalb von chart)
virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point<3> & /*p3*/) const
{ gi.trignum = 1; return 1;};
virtual int CalcPointGeomInfo(int /* surfind */, PointGeomInfo& gi, const Point<3> & p3) const
{ return CalcPointGeomInfo (gi, p3); }
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & gi, Vec<3> & n) const;
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const;
void CheckMeshApproximation (Mesh & mesh);

View File

@ -19,15 +19,16 @@ namespace netgen
};
void MeshOptimize2d :: GenericImprove (Mesh & mesh)
void MeshOptimize2d :: GenericImprove ()
{
static Timer timer("MeshOptimize2d::GenericImprove"); RegionTimer reg(timer);
if (!faceindex)
{
if (writestatus)
PrintMessage (3, "Generic Improve");
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
GenericImprove (mesh);
GenericImprove ();
faceindex = 0;
}
@ -395,10 +396,8 @@ namespace netgen
// calc metric badness
double bad1 = 0, bad2 = 0;
Vec<3> n;
SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1));
GetNormalVector (surfnr, mesh.Point(pmap.Get(1)), pgi.Elem(1), n);
// SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1));
auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), &pgi.Elem(1));
for (int j = 0; j < rule.oldels.Size(); j++)
bad1 += mesh[elmap[j]].CalcJacobianBadness (mesh.Points(), n);

View File

@ -183,7 +183,7 @@ void MeshOptimize3d :: CombineImproveSequential (Mesh & mesh,
// mesh.CalcSurfacesOfNode ();
const char * savetask = multithread.task;
multithread.task = "Combine Improve";
multithread.task = "Optimize Volume: Combine Improve";
double totalbad = 0;
@ -198,7 +198,7 @@ void MeshOptimize3d :: CombineImproveSequential (Mesh & mesh,
if (goal == OPT_QUALITY)
{
totalbad = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
totalbad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << totalbad << endl;
PrintMessage (5, "Total badness = ", totalbad);
}
@ -395,7 +395,7 @@ void MeshOptimize3d :: CombineImproveSequential (Mesh & mesh,
if (goal == OPT_QUALITY)
{
totalbad = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
totalbad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << totalbad << endl;
int cntill = 0;
@ -409,60 +409,6 @@ void MeshOptimize3d :: CombineImproveSequential (Mesh & mesh,
multithread.task = savetask;
}
void MeshOptimize3d :: BuildEdgeList( const Mesh & mesh, const Table<ElementIndex, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges )
{
static Timer tbuild_edges("Build edges"); RegionTimer reg(tbuild_edges);
static constexpr int tetedges[6][2] =
{ { 0, 1 }, { 0, 2 }, { 0, 3 },
{ 1, 2 }, { 1, 3 }, { 2, 3 } };
int ntasks = 2*ngcore::TaskManager::GetMaxThreads();
Array<Array<std::tuple<PointIndex,PointIndex>>> task_edges(ntasks);
ParallelFor(IntRange(ntasks), [&] (int ti)
{
auto myrange = mesh.Points().Range().Split(ti, ntasks);
ArrayMem<std::tuple<PointIndex,PointIndex>, 100> local_edges;
for (auto pi : myrange)
{
local_edges.SetSize(0);
for(auto ei : elementsonnode[pi])
{
const Element & elem = mesh[ei];
if (elem.IsDeleted()) continue;
for (int j = 0; j < 6; j++)
{
PointIndex pi0 = elem[tetedges[j][0]];
PointIndex pi1 = elem[tetedges[j][1]];
if (pi1 < pi0) Swap(pi0, pi1);
if(pi0==pi)
local_edges.Append(std::make_tuple(pi0, pi1));
}
}
QuickSort(local_edges);
auto edge_prev = std::make_tuple<PointIndex, PointIndex>(-1,-1);
for(auto edge : local_edges)
if(edge != edge_prev)
{
task_edges[ti].Append(edge);
edge_prev = edge;
}
}
}, ntasks);
int num_edges = 0;
for (auto & edg : task_edges)
num_edges += edg.Size();
edges.SetAllocSize(num_edges);
for (auto & edg : task_edges)
edges.Append(edg);
}
void MeshOptimize3d :: CombineImprove (Mesh & mesh,
OPTIMIZEGOAL goal)
{
@ -489,7 +435,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
// mesh.CalcSurfacesOfNode ();
const char * savetask = multithread.task;
multithread.task = "Combine Improve";
multithread.task = "Optimize Volume: Combine Improve";
tbad.Start();
@ -511,7 +457,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
if (goal == OPT_QUALITY)
{
totalbad = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
totalbad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << totalbad << endl;
}
@ -565,7 +511,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
if (goal == OPT_QUALITY)
{
totalbad = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
totalbad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << totalbad << endl;
int cntill = 0;
@ -581,6 +527,275 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table<ElementIndex,PointIndex> & elementsonnode, Array<double> &elerrs, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only)
{
double d_badness = 0.0;
int cnt = 0;
ArrayMem<ElementIndex, 20> hasbothpoints;
if (mesh.BoundaryEdge (pi1, pi2)) return 0.0;
for (ElementIndex ei : elementsonnode[pi1])
{
Element & el = mesh[ei];
if(el.IsDeleted()) return 0.0;
if (mesh[ei].GetType() != TET) return 0.0;
bool has1 = el.PNums().Contains(pi1);
bool has2 = el.PNums().Contains(pi2);
if (has1 && has2)
if (!hasbothpoints.Contains (ei))
hasbothpoints.Append (ei);
}
if(mp.only3D_domain_nr)
for(auto ei : hasbothpoints)
if(mp.only3D_domain_nr != mesh[ei].GetIndex())
return 0.0;
if (goal == OPT_LEGAL)
{
bool all_tets_legal = true;
for(auto ei : hasbothpoints)
if( !mesh.LegalTet (mesh[ei]) || elerrs[ei] > 1e3)
all_tets_legal = false;
if(all_tets_legal)
return 0.0;
}
double bad1 = 0.0;
double bad1_max = 0.0;
for (ElementIndex ei : hasbothpoints)
{
double bad = elerrs[ei];
bad1 += bad;
bad1_max = max(bad1_max, bad);
}
if(bad1_max < 100.0)
return 0.0;
bool puretet = 1;
for (ElementIndex ei : hasbothpoints)
if (mesh[ei].GetType() != TET)
puretet = 0;
if (!puretet) return 0.0;
Point3d p1 = mesh[pi1];
Point3d p2 = mesh[pi2];
locfaces.SetSize(0);
for (ElementIndex ei : hasbothpoints)
{
const Element & el = mesh[ei];
for (int l = 0; l < 4; l++)
if (el[l] == pi1 || el[l] == pi2)
{
INDEX_3 i3;
Element2d face(TRIG);
el.GetFace (l+1, face);
for (int kk = 1; kk <= 3; kk++)
i3.I(kk) = face.PNum(kk);
locfaces.Append (i3);
}
}
PointFunction1 pf (mesh.Points(), locfaces, mp, -1);
OptiParameters par;
par.maxit_linsearch = 50;
par.maxit_bfgs = 20;
Point3d pnew = Center (p1, p2);
Vector px(3);
px(0) = pnew.X();
px(1) = pnew.Y();
px(2) = pnew.Z();
if (bad1_max > 0.1 * badmax)
{
int pok = pf.Func (px) < 1e10;
if (!pok)
pok = FindInnerPoint (mesh.Points(), locfaces, pnew);
if(pok)
{
px(0) = pnew.X();
px(1) = pnew.Y();
px(2) = pnew.Z();
BFGS (px, pf, par);
pnew.X() = px(0);
pnew.Y() = px(1);
pnew.Z() = px(2);
}
}
double bad2 = pf.Func (px);
mesh[ptmp] = Point<3>(pnew);
for (int k = 0; k < hasbothpoints.Size(); k++)
{
Element & oldel = mesh[hasbothpoints[k]];
Element newel1 = oldel;
Element newel2 = oldel;
oldel.flags.illegal_valid = 0;
newel1.flags.illegal_valid = 0;
newel2.flags.illegal_valid = 0;
for (int l = 0; l < 4; l++)
{
if (newel1[l] == pi2) newel1[l] = ptmp;
if (newel2[l] == pi1) newel2[l] = ptmp;
}
if (!mesh.LegalTet (oldel)) bad1 += 1e6;
if (!mesh.LegalTet (newel1)) bad2 += 1e6;
if (!mesh.LegalTet (newel2)) bad2 += 1e6;
}
d_badness = bad2-bad1;
if(check_only)
return d_badness;
if (d_badness<0.0)
{
cnt++;
PointIndex pinew = mesh.AddPoint (pnew);
for (ElementIndex ei : hasbothpoints)
{
Element & oldel = mesh[ei];
Element newel1 = oldel;
Element newel2 = oldel;
oldel.flags.illegal_valid = 0;
oldel.Delete();
newel1.flags.illegal_valid = 0;
newel2.flags.illegal_valid = 0;
for (int l = 0; l < 4; l++)
{
if (newel1[l] == pi2) newel1[l] = pinew;
if (newel2[l] == pi1) newel2[l] = pinew;
}
mesh.AddVolumeElement (newel1);
mesh.AddVolumeElement (newel2);
}
}
return d_badness;
}
void MeshOptimize3d :: SplitImprove (Mesh & mesh,
OPTIMIZEGOAL goal)
{
static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t);
static Timer topt("Optimize");
static Timer tsearch("Search");
// return SplitImproveSequential(mesh, goal);
int np = mesh.GetNP();
int ne = mesh.GetNE();
double bad = 0.0;
double badmax = 0.0;
auto elementsonnode = mesh.CreatePoint2ElementTable();
Array<double> elerrs(ne);
const char * savetask = multithread.task;
multithread.task = "Optimize Volume: Split Improve";
PrintMessage (3, "SplitImprove");
(*testout) << "start SplitImprove" << "\n";
mesh.BoundaryEdge (1,2); // ensure the boundary-elements table is built
ParallelFor( mesh.VolumeElements().Range(), [&] (ElementIndex ei) NETGEN_LAMBDA_INLINE
{
if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
return;
elerrs[ei] = CalcBad (mesh.Points(), mesh[ei], 0);
bad += elerrs[ei];
AtomicMax(badmax, elerrs[ei]);
});
if (goal == OPT_QUALITY)
{
bad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad << endl;
}
Array<std::tuple<PointIndex,PointIndex>> edges;
BuildEdgeList(mesh, elementsonnode, edges);
// Find edges with improvement
Array<std::tuple<double, int>> candidate_edges(edges.Size());
std::atomic<int> improvement_counter(0);
auto ptmp = mesh.AddPoint( {0,0,0} );
tsearch.Start();
ParallelForRange(Range(edges), [&] (auto myrange)
{
NgArray<INDEX_3> locfaces;
for(auto i : myrange)
{
auto [p0,p1] = edges[i];
double d_badness = SplitImproveEdge (mesh, goal, elementsonnode, elerrs, locfaces, badmax, p0, p1, ptmp, true);
if(d_badness<0.0)
{
int index = improvement_counter++;
candidate_edges[index] = make_tuple(d_badness, i);
}
}
}, ngcore::TasksPerThread(4));
tsearch.Stop();
auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
QuickSort(edges_with_improvement);
PrintMessage(5, edges.Size(), " edges");
PrintMessage(5, edges_with_improvement.Size(), " edges with improvement");
// Apply actual optimizations
topt.Start();
int cnt = 0;
NgArray<INDEX_3> locfaces;
for(auto [d_badness, ei] : edges_with_improvement)
{
auto [p0,p1] = edges[ei];
if (SplitImproveEdge (mesh, goal, elementsonnode, elerrs, locfaces, badmax, p0, p1, ptmp, false) < 0.0)
cnt++;
}
topt.Stop();
mesh.Compress();
PrintMessage (5, cnt, " splits performed");
(*testout) << "Splitt - Improve done" << "\n";
if (goal == OPT_QUALITY)
{
bad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad << endl;
int cntill = 0;
ne = mesh.GetNE();
for (ElementIndex ei = 0; ei < ne; ei++)
if (!mesh.LegalTet (mesh[ei]))
cntill++;
// cout << cntill << " illegal tets" << endl;
}
multithread.task = savetask;
}
/*
@ -588,7 +803,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
If mesh quality is improved by inserting a node into an inner edge,
the edge is split into two parts.
*/
void MeshOptimize3d :: SplitImprove (Mesh & mesh,
void MeshOptimize3d :: SplitImproveSequential (Mesh & mesh,
OPTIMIZEGOAL goal)
{
static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t);
@ -611,7 +826,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
illegaltet.Clear();
const char * savetask = multithread.task;
multithread.task = "Split Improve";
multithread.task = "Optimize Volume: Split Improve";
PrintMessage (3, "SplitImprove");
(*testout) << "start SplitImprove" << "\n";
@ -642,7 +857,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
if (goal == OPT_QUALITY)
{
bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
}
@ -864,7 +1079,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
if (goal == OPT_QUALITY)
{
bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
int cntill = 0;
@ -906,7 +1121,7 @@ void MeshOptimize3d :: SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal,
(*testout) << "\n" << "Start SwapImprove" << endl;
const char * savetask = multithread.task;
multithread.task = "Swap Improve";
multithread.task = "Optimize Volume: Swap Improve";
// mesh.CalcSurfacesOfNode ();
@ -925,7 +1140,7 @@ void MeshOptimize3d :: SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal,
// Calculate total badness
if (goal == OPT_QUALITY)
{
bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
}
@ -1752,7 +1967,7 @@ void MeshOptimize3d :: SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal,
}
bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
const NgBitArray * working_elements,
Table<ElementIndex, PointIndex> & elementsonnode,
INDEX_3_HASHTABLE<int> & faces,
@ -1768,10 +1983,10 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
Element el1b(TET), el2b(TET), el3b(TET), el4b(TET);
ArrayMem<ElementIndex, 20> hasbothpoints;
bool do_swap = false;
double d_badness = 0.0;
if (pi2 < pi1) Swap (pi1, pi2);
if (mesh.BoundaryEdge (pi1, pi2)) return false;
if (mesh.BoundaryEdge (pi1, pi2)) return 0.0;
hasbothpoints.SetSize (0);
@ -1780,7 +1995,7 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
bool has1 = 0, has2 = 0;
const Element & elem = mesh[elnr];
if (elem.IsDeleted()) return false;
if (elem.IsDeleted()) return 0.0;
for (int l = 0; l < elem.GetNP(); l++)
{
@ -1803,27 +2018,27 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
for (ElementIndex ei : hasbothpoints)
{
if (mesh[ei].GetType () != TET)
return false;
return 0.0;
if (mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
return false;
return 0.0;
if ((mesh.ElementType(ei)) == FIXEDELEMENT)
return false;
return 0.0;
if(working_elements &&
ei < working_elements->Size() &&
!working_elements->Test(ei))
return false;
return 0.0;
if (mesh[ei].IsDeleted())
return false;
return 0.0;
if ((goal == OPT_LEGAL) &&
mesh.LegalTet (mesh[ei]) &&
CalcBad (mesh.Points(), mesh[ei], 0) < 1e3)
return false;
return 0.0;
}
int nsuround = hasbothpoints.Size();
@ -1937,8 +2152,9 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
{
// (*mycout) << "3->2 " << flush;
// (*testout) << "3->2 conversion" << endl;
do_swap = true;
if(check_only) return do_swap;
d_badness = bad2-bad1;
if(check_only)
return d_badness;
/*
@ -2133,12 +2349,9 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
swap3 = !swap2 && (bad3 < bad1);
}
if (swap2 || swap3)
{
do_swap = true;
if(check_only) return do_swap;
}
d_badness = swap2 ? bad2-bad1 : bad3-bad1;
if(check_only)
return d_badness;
if (swap2)
{
@ -2333,8 +2546,9 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
if (bestl != -1)
{
// (*mycout) << nsuround << "->" << 2 * (nsuround-2) << " " << flush;
do_swap = true;
if(check_only) return do_swap;
d_badness = badopt-bad1;
if(check_only)
return d_badness;
for (int k = bestl+1; k <= nsuround + bestl - 2; k++)
{
@ -2377,7 +2591,7 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
}
}
}
return do_swap;
return d_badness;
}
void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
@ -2403,7 +2617,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
(*testout) << "\n" << "Start SwapImprove" << endl;
const char * savetask = multithread.task;
multithread.task = "Swap Improve";
multithread.task = "Optimize Volume: Swap Improve";
INDEX_3_HASHTABLE<int> faces(mesh.GetNOpenElements()/3 + 2);
if (goal == OPT_CONFORM)
@ -2420,14 +2634,14 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
// Calculate total badness
if (goal == OPT_QUALITY)
{
double bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
double bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
}
Array<std::tuple<PointIndex,PointIndex>> edges;
BuildEdgeList(mesh, elementsonnode, edges);
Array<int> candidate_edges(edges.Size());
Array<std::tuple<double, int>> candidate_edges(edges.Size());
std::atomic<int> improvement_counter(0);
tloop.Start();
@ -2440,18 +2654,22 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
break;
auto [pi0, pi1] = edges[i];
if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, true))
candidate_edges[improvement_counter++] = i;
double d_badness = SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, true);
if(d_badness<0.0)
{
int index = improvement_counter++;
candidate_edges[index] = make_tuple(d_badness, i);
}
}
});
auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
QuickSort(edges_with_improvement);
for(auto ei : edges_with_improvement)
for(auto [d_badness, ei] : edges_with_improvement)
{
auto [pi0,pi1] = edges[ei];
if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, false))
if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, false) < 0.0)
cnt++;
}
@ -3331,7 +3549,7 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
2 -> 3 conversion
*/
bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face,
double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face,
Table<ElementIndex, PointIndex> & elementsonnode,
TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only )
{
@ -3339,9 +3557,10 @@ bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementInd
Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
int j = face;
double bad1, bad2;
double d_badness = 0.0;
Element & elem = mesh[eli1];
if (elem.IsDeleted()) return false;
if (elem.IsDeleted()) return 0.0;
int mattyp = elem.GetIndex();
@ -3387,7 +3606,7 @@ bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementInd
}
}
if (bface) return false;
if (bface) return 0.0;
FlatArray<ElementIndex> row = elementsonnode[pi1];
@ -3457,14 +3676,16 @@ bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementInd
bad2 += 1e4;
bool do_swap = (bad2 < bad1);
d_badness = bad2 - bad1;
if ( ((bad2 < 1e6) || (bad2 < 10 * bad1)) &&
mesh.BoundaryEdge (pi4, pi5))
do_swap = 1;
d_badness = -1e4;
if(check_only)
return d_badness;
if (!check_only && do_swap)
if (d_badness<0.0)
{
el31.flags.illegal_valid = 0;
el32.flags.illegal_valid = 0;
@ -3476,12 +3697,11 @@ bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementInd
mesh.AddVolumeElement (el32);
mesh.AddVolumeElement (el33);
}
return do_swap;
return d_badness;
}
}
}
return false;
return d_badness;
}
/*
@ -3527,7 +3747,7 @@ void MeshOptimize3d :: SwapImprove2Sequential (Mesh & mesh, OPTIMIZEGOAL goal)
// Calculate total badness
bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
cout << "tot bad = " << bad1 << endl;
@ -3582,7 +3802,8 @@ void MeshOptimize3d :: SwapImprove2Sequential (Mesh & mesh, OPTIMIZEGOAL goal)
*/
bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
mesh.Compress();
bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
(*testout) << "swapimprove2 done" << "\n";
// (*mycout) << "Vol = " << CalcVolume (points, volelements) << "\n";
@ -3611,7 +3832,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
PrintMessage (3, "SwapImprove2 ");
(*testout) << "\n" << "Start SwapImprove2" << "\n";
bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
// find elements on node
@ -3624,8 +3845,8 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
int num_threads = ngcore::TaskManager::GetNumThreads();
Array<std::tuple<ElementIndex, int>> faces_with_improvement;
Array<Array<std::tuple<ElementIndex, int>>> faces_with_improvement_threadlocal(num_threads);
Array<std::tuple<double, ElementIndex, int>> faces_with_improvement;
Array<Array<std::tuple<double, ElementIndex, int>>> faces_with_improvement_threadlocal(num_threads);
ParallelForRange( Range(ne), [&]( auto myrange )
{
@ -3652,8 +3873,9 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
for (int j = 0; j < 4; j++)
{
if(SwapImprove2( mesh, goal, eli1, j, elementsonnode, belementsonnode, true))
my_faces_with_improvement.Append( std::make_tuple(eli1, j) );
double d_badness = SwapImprove2( mesh, goal, eli1, j, elementsonnode, belementsonnode, true);
if(d_badness<0.0)
my_faces_with_improvement.Append( std::make_tuple(d_badness, eli1, j) );
}
}
});
@ -3663,12 +3885,14 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
QuickSort(faces_with_improvement);
for (auto [eli,j] : faces_with_improvement)
cnt += SwapImprove2( mesh, goal, eli, j, elementsonnode, belementsonnode, false);
for (auto [dummy, eli,j] : faces_with_improvement)
if(SwapImprove2( mesh, goal, eli, j, elementsonnode, belementsonnode, false) < 0.0)
cnt++;
PrintMessage (5, cnt, " swaps performed");
bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
mesh.Compress();
bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
(*testout) << "swapimprove2 done" << "\n";
}

View File

@ -12,8 +12,6 @@ class MeshOptimize3d
{
const MeshingParameters & mp;
void BuildEdgeList( const Mesh & mesh, const Table<ElementIndex, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges );
public:
MeshOptimize3d (const MeshingParameters & amp) : mp(amp) { ; }
@ -26,8 +24,11 @@ public:
void CombineImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
void SplitImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
double SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table<ElementIndex,PointIndex> & elementsonnode, Array<double> &elerrs, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false);
bool SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, const NgBitArray * working_elements, Table<ElementIndex,PointIndex> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false);
double SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, const NgBitArray * working_elements, Table<ElementIndex,PointIndex> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false);
void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
const NgBitArray * working_elements = NULL);
void SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
@ -37,7 +38,7 @@ public:
const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL);
void SwapImprove2Sequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
void SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
bool SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only=false );
double SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only=false );
double
CalcBad (const Mesh::T_POINTS & points, const Element & elem, double h)

View File

@ -1,7 +1,11 @@
#include <mystdlib.h>
#include <atomic>
#include "meshing.hpp"
#ifdef NG_PYTHON
// must be included to instantiate Archive::Shallow(NetgenGeometry&)
#include <core/python_ngcore.hpp>
#endif
namespace netgen
{
@ -1662,6 +1666,8 @@ namespace netgen
void Mesh :: CalcSurfacesOfNode ()
{
static Timer t("Mesh::CalcSurfacesOfNode"); RegionTimer reg (t);
static Timer tn2se("Mesh::CalcSurfacesOfNode - surf on node");
static Timer tht("Mesh::CalcSurfacesOfNode - surfelementht");
// surfacesonnode.SetSize (GetNP());
TABLE<int,PointIndex::BASE> surfacesonnode(GetNP());
@ -1683,6 +1689,7 @@ namespace netgen
surfelementht = make_unique<INDEX_3_CLOSED_HASHTABLE<int>> (3*GetNSE() + 1);
segmentht = make_unique<INDEX_2_CLOSED_HASHTABLE<int>> (3*GetNSeg() + 1);
tn2se.Start();
if (dimension == 3)
/*
for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
@ -1734,7 +1741,9 @@ namespace netgen
surfelementht -> AllocateElements();
*/
tn2se.Stop();
tht.Start();
if (dimension==3)
for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
{
@ -1748,6 +1757,7 @@ namespace netgen
i3.Sort();
surfelementht -> Set (i3, sei); // war das wichtig ??? sel.GetIndex());
}
tht.Stop();
// int np = GetNP();
@ -1878,10 +1888,43 @@ namespace netgen
int ne = GetNE();
int nse = GetNSE();
NgArray<int,PointIndex::BASE> numonpoint(np);
numonpoint = 0;
t_table.Start();
TableCreator<ElementIndex, PointIndex> creator(np);
for ( ; !creator.Done(); creator++)
// for (ElementIndex ei : Range(VolumeElements()))
ParallelFor
(Range(VolumeElements()), [&] (ElementIndex ei)
{
const Element & el = (*this)[ei];
if (dom == 0 || dom == el.GetIndex())
{
if (el.GetNP() == 4)
{
INDEX_4 i4(el[0], el[1], el[2], el[3]);
i4.Sort();
creator.Add (PointIndex(i4.I1()), ei);
creator.Add (PointIndex(i4.I2()), ei);
}
else
{
for (PointIndex pi : el.PNums())
creator.Add(pi, ei);
}
}
});
auto elsonpoint = creator.MoveTable();
ParallelFor (Range(elsonpoint), [&] (auto i)
{
QuickSort(elsonpoint[i]);
});
NgArray<int,PointIndex::BASE> numonpoint(np);
/*
numonpoint = 0;
for (ElementIndex ei = 0; ei < ne; ei++)
{
const Element & el = (*this)[ei];
@ -1918,6 +1961,7 @@ namespace netgen
elsonpoint.Add (el[j], ei);
}
}
*/
t_table.Stop();
@ -3688,63 +3732,50 @@ namespace netgen
{
static Timer t("Mesh::CheckOverlappingBoundary"); RegionTimer reg(t);
int i, j, k;
Point3d pmin, pmax;
GetBox (pmin, pmax);
BoxTree<3> setree(pmin, pmax);
NgArray<int> inters;
BoxTree<3, SurfaceElementIndex> setree(pmin, pmax);
// NgArray<SurfaceElementIndex> inters;
bool overlap = 0;
bool incons_layers = 0;
/*
for (i = 1; i <= GetNSE(); i++)
SurfaceElement(i).badel = 0;
*/
for (Element2d & el : SurfaceElements())
el.badel = false;
for (i = 1; i <= GetNSE(); i++)
for (SurfaceElementIndex sei : Range(SurfaceElements()))
{
const Element2d & tri = SurfaceElement(i);
const Element2d & tri = SurfaceElement(sei);
Point3d tpmin (Point(tri[0]));
Point3d tpmax (tpmin);
Box<3> box(Box<3>::EMPTY_BOX);
for (PointIndex pi : tri.PNums())
box.Add (Point(pi));
for (k = 1; k < tri.GetNP(); k++)
{
tpmin.SetToMin (Point (tri[k]));
tpmax.SetToMax (Point (tri[k]));
}
Vec3d diag(tpmin, tpmax);
tpmax = tpmax + 0.1 * diag;
tpmin = tpmin - 0.1 * diag;
setree.Insert (tpmin, tpmax, i);
box.Increase(1e-3*box.Diam());
setree.Insert (box, sei);
}
for (i = 1; i <= GetNSE(); i++)
std::mutex m;
// for (SurfaceElementIndex sei : Range(SurfaceElements()))
ParallelForRange
(Range(SurfaceElements()), [&] (auto myrange)
{
const Element2d & tri = SurfaceElement(i);
Point3d tpmin (Point(tri[0]));
Point3d tpmax (tpmin);
for (k = 1; k < tri.GetNP(); k++)
for (SurfaceElementIndex sei : myrange)
{
tpmin.SetToMin (Point (tri[k]));
tpmax.SetToMax (Point (tri[k]));
}
const Element2d & tri = SurfaceElement(sei);
setree.GetIntersecting (tpmin, tpmax, inters);
Box<3> box(Box<3>::EMPTY_BOX);
for (PointIndex pi : tri.PNums())
box.Add (Point(pi));
for (j = 1; j <= inters.Size(); j++)
setree.GetFirstIntersecting
(box.PMin(), box.PMax(),
[&] (SurfaceElementIndex sej)
{
const Element2d & tri2 = SurfaceElement(inters.Get(j));
const Element2d & tri2 = SurfaceElement(sej);
if ( (*this)[tri[0]].GetLayer() != (*this)[tri2[0]].GetLayer())
continue;
return false;
if ( (*this)[tri[0]].GetLayer() != (*this)[tri[1]].GetLayer() ||
(*this)[tri[0]].GetLayer() != (*this)[tri[2]].GetLayer())
@ -3753,66 +3784,52 @@ namespace netgen
cout << "inconsistent layers in triangle" << endl;
}
const netgen::Point<3> *trip1[3], *trip2[3];
for (k = 1; k <= 3; k++)
for (int k = 0; k < 3; k++)
{
trip1[k-1] = &Point (tri.PNum(k));
trip2[k-1] = &Point (tri2.PNum(k));
trip1[k] = &Point (tri[k]);
trip2[k] = &Point (tri2[k]);
}
if (IntersectTriangleTriangle (&trip1[0], &trip2[0]))
{
overlap = 1;
lock_guard<std::mutex> guard(m);
PrintWarning ("Intersecting elements "
,i, " and ", inters.Get(j));
,int(sei), " and ", int(sej));
(*testout) << "Intersecting: " << endl;
(*testout) << "openelement " << i << " with open element " << inters.Get(j) << endl;
(*testout) << "openelement " << sei << " with open element " << sej << endl;
cout << "el1 = " << tri << endl;
cout << "el2 = " << tri2 << endl;
cout << "layer1 = " << (*this)[tri[0]].GetLayer() << endl;
cout << "layer2 = " << (*this)[tri2[0]].GetLayer() << endl;
for (k = 1; k <= 3; k++)
for (int k = 1; k <= 3; k++)
(*testout) << tri.PNum(k) << " ";
(*testout) << endl;
for (k = 1; k <= 3; k++)
for (int k = 1; k <= 3; k++)
(*testout) << tri2.PNum(k) << " ";
(*testout) << endl;
for (k = 0; k <= 2; k++)
for (int k = 0; k <= 2; k++)
(*testout) << *trip1[k] << " ";
(*testout) << endl;
for (k = 0; k <= 2; k++)
for (int k = 0; k <= 2; k++)
(*testout) << *trip2[k] << " ";
(*testout) << endl;
(*testout) << "Face1 = " << GetFaceDescriptor(tri.GetIndex()) << endl;
(*testout) << "Face1 = " << GetFaceDescriptor(tri2.GetIndex()) << endl;
/*
INDEX_3 i3(tri.PNum(1), tri.PNum(2), tri.PNum(3));
i3.Sort();
for (k = 1; k <= GetNSE(); k++)
{
const Element2d & el2 = SurfaceElement(k);
INDEX_3 i3b(el2.PNum(1), el2.PNum(2), el2.PNum(3));
i3b.Sort();
if (i3 == i3b)
{
SurfaceElement(k).badel = 1;
SurfaceElement(sei).badel = 1;
SurfaceElement(sej).badel = 1;
}
return false;
});
}
*/
SurfaceElement(i).badel = 1;
SurfaceElement(inters.Get(j)).badel = 1;
}
}
}
});
// bug 'fix'
if (incons_layers) overlap = 0;
@ -3924,6 +3941,44 @@ namespace netgen
return 1;
}
double Mesh :: CalcTotalBad (const MeshingParameters & mp )
{
static Timer t("CalcTotalBad"); RegionTimer reg(t);
static constexpr int n_classes = 20;
double sum = 0;
tets_in_qualclass.SetSize(n_classes);
tets_in_qualclass = 0;
ParallelForRange( IntRange(volelements.Size()), [&] (auto myrange)
{
double local_sum = 0.0;
double teterrpow = mp.opterrpow;
std::array<int,n_classes> classes_local{};
for (auto i : myrange)
{
double elbad = pow (max2(CalcBad (points, volelements[i], 0, mp),1e-10), 1/teterrpow);
int qualclass = int (n_classes / elbad + 1);
if (qualclass < 1) qualclass = 1;
if (qualclass > n_classes) qualclass = n_classes;
classes_local[qualclass-1]++;
local_sum += elbad;
}
AtomicAdd(sum, local_sum);
for (auto i : Range(n_classes))
AsAtomic(tets_in_qualclass[i]) += classes_local[i];
});
return sum;
}
@ -6154,14 +6209,19 @@ namespace netgen
{
for (PointIndex pi : myrange)
QuickSort(elementsonnode[pi]);
});
}, ngcore::TasksPerThread(4));
return move(elementsonnode);
}
Table<SurfaceElementIndex, PointIndex> Mesh :: CreatePoint2SurfaceElementTable() const
Table<SurfaceElementIndex, PointIndex> Mesh :: CreatePoint2SurfaceElementTable( int faceindex ) const
{
static Timer timer("Mesh::CreatePoint2SurfaceElementTable"); RegionTimer rt(timer);
TableCreator<SurfaceElementIndex, PointIndex> creator(GetNP());
if(faceindex==0)
{
for ( ; !creator.Done(); creator++)
ngcore::ParallelForRange
(Range(surfelements), [&] (auto myrange)
@ -6169,7 +6229,21 @@ namespace netgen
for (SurfaceElementIndex ei : myrange)
for (PointIndex pi : (*this)[ei].PNums())
creator.Add (pi, ei);
});
}, ngcore::TasksPerThread(4));
}
else
{
Array<SurfaceElementIndex> face_els;
GetSurfaceElementsOfFace(faceindex, face_els);
for ( ; !creator.Done(); creator++)
ngcore::ParallelForRange
(Range(face_els), [&] (auto myrange)
{
for (auto i : myrange)
for (PointIndex pi : (*this)[face_els[i]].PNums())
creator.Add (pi, face_els[i]);
}, ngcore::TasksPerThread(4));
}
auto elementsonnode = creator.MoveTable();
ngcore::ParallelForRange
@ -6512,6 +6586,21 @@ namespace netgen
return defaultstring;
}
NgArray<string*> & Mesh :: GetRegionNamesCD (int codim)
{
switch (codim)
{
case 0: return materials;
case 1: return bcnames;
case 2: return cd2names;
case 3: return cd3names;
default: throw Exception("don't have regions of co-dimension "+ToString(codim));
}
}
void Mesh :: SetUserData(const char * id, NgArray<int> & data)
{
if(userdata_int.Used(id))

View File

@ -62,6 +62,8 @@ namespace netgen
/// open segmenets for surface meshing
NgArray<Segment> opensegments;
Array<int> tets_in_qualclass;
/**
@ -559,6 +561,10 @@ namespace netgen
*/
void FreeOpenElementsEnvironment (int layers);
DLL_HEADER double CalcTotalBad (const MeshingParameters & mp);
FlatArray<int> GetQualityHistogram() { return tets_in_qualclass; }
///
bool LegalTet (Element & el) const
{
@ -681,6 +687,9 @@ namespace netgen
string * GetBCNamePtr (int bcnr) const
{ return (bcnr < bcnames.Size() && bcnames[bcnr]) ? bcnames[bcnr] : &default_bc; }
NgArray<string*> & GetRegionNamesCD (int codim);
///
void ClearFaceDescriptors()
{ facedecoding.SetSize(0); }
@ -753,7 +762,7 @@ namespace netgen
Table<ElementIndex, PointIndex> CreatePoint2ElementTable() const;
Table<SurfaceElementIndex, PointIndex> CreatePoint2SurfaceElementTable() const;
Table<SurfaceElementIndex, PointIndex> CreatePoint2SurfaceElementTable( int faceindex=0 ) const;
DLL_HEADER bool PureTrigMesh (int faceindex = 0) const;
DLL_HEADER bool PureTetMesh () const;
@ -855,7 +864,7 @@ namespace netgen
///
friend class Meshing3;
// only for saving the geometry
enum GEOM_TYPE { NO_GEOM = 0, GEOM_2D = 1, GEOM_CSG = 10, GEOM_STL = 11, GEOM_OCC = 12, GEOM_ACIS = 13 };
GEOM_TYPE geomtype;

View File

@ -646,6 +646,8 @@ namespace netgen
{
static Timer t("OptimizeVolume"); RegionTimer reg(t);
RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0);
const char* savetask = multithread.task;
multithread.task = "Optimize Volume";
int i;
@ -663,7 +665,7 @@ namespace netgen
*/
mesh3d.CalcSurfacesOfNode();
for (i = 1; i <= mp.optsteps3d; i++)
for (auto i : Range(mp.optsteps3d))
{
if (multithread.terminate)
break;
@ -672,12 +674,13 @@ namespace netgen
// teterrpow = mp.opterrpow;
// for (size_t j = 1; j <= strlen(mp.optimize3d); j++)
for (size_t j = 1; j <= mp.optimize3d.length(); j++)
for (auto j : Range(mp.optimize3d.size()))
{
multithread.percent = 100.* (double(j)/mp.optimize3d.size() + i)/mp.optsteps3d;
if (multithread.terminate)
break;
switch (mp.optimize3d[j-1])
switch (mp.optimize3d[j])
{
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
case 'd': optmesh.SplitImprove(mesh3d); break;
@ -698,6 +701,7 @@ namespace netgen
MeshQuality3d (mesh3d);
}
multithread.task = savetask;
return MESHING3_OK;
}

View File

@ -31,30 +31,30 @@ namespace netgen
{
case 's':
{ // topological swap
MeshOptimize2d meshopt;
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.EdgeSwapping (mesh, 0);
meshopt.EdgeSwapping (0);
break;
}
case 'S':
{ // metric swap
MeshOptimize2d meshopt;
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.EdgeSwapping (mesh, 1);
meshopt.EdgeSwapping (1);
break;
}
case 'm':
{
MeshOptimize2d meshopt;
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.ImproveMesh(mesh, mp);
meshopt.ImproveMesh(mp);
break;
}
case 'c':
{
MeshOptimize2d meshopt;
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.CombineImprove(mesh);
meshopt.CombineImprove();
break;
}
default:

View File

@ -42,13 +42,12 @@ namespace netgen
#define _INCLUDE_MORE
#include "findip.hpp"
#include "findip2.hpp"
#include "meshing3.hpp"
#include "improve3.hpp"
#include "findip.hpp"
#include "findip2.hpp"
#include "curvedelems.hpp"
#include "clusters.hpp"

View File

@ -6,6 +6,14 @@
namespace netgen
{
ostream& operator << (ostream& ost, const MultiPointGeomInfo& mpgi)
{
for(auto i : Range(mpgi.GetNPGI()))
{
ost << "gi[" << i << "] = " << mpgi.GetPGI(i+1) << endl;
}
return ost;
}
static void glrender (int wait);
#ifdef OPENGL
extern DLL_HEADER void Render(bool blocking = false);
@ -30,8 +38,10 @@ namespace netgen
static Array<unique_ptr<netrule>> global_quad_rules;
Meshing2 :: Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox)
: adfront(aboundingbox), boundingbox(aboundingbox)
Meshing2 :: Meshing2 (const NetgenGeometry& ageo,
const MeshingParameters & mp,
const Box<3> & aboundingbox)
: geo(ageo), adfront(aboundingbox), boundingbox(aboundingbox)
{
static Timer t("Mesing2::Meshing2"); RegionTimer r(t);
@ -125,28 +135,37 @@ namespace netgen
// static Vec3d ex, ey;
// static Point3d globp1;
void Meshing2 :: DefineTransformation (const Point<3> & p1, const Point<3> & p2,
const PointGeomInfo * geominfo1,
const PointGeomInfo * geominfo2)
void Meshing2 :: DefineTransformation (const Point<3> & ap1,
const Point<3> & ap2,
const PointGeomInfo * gi1,
const PointGeomInfo * gi2)
{
globp1 = p1;
ex = p2 - p1;
ex /= ex.Length();
ey.X() = -ex.Y();
ey.Y() = ex.X();
ey.Z() = 0;
p1 = ap1;
p2 = ap2;
auto n1 = geo.GetNormal(gi1->trignum, p1, gi1);
auto n2 = geo.GetNormal(gi2->trignum, p2, gi2);
ez = 0.5 * (n1+n2);
ez.Normalize();
ex = (p2-p1).Normalize();
ez -= (ez*ex)*ex;
ez.Normalize();
ey = Cross(ez, ex);
}
void Meshing2 :: TransformToPlain (const Point<3> & locpoint,
const MultiPointGeomInfo & geominf,
const MultiPointGeomInfo & geominfo,
Point<2> & plainpoint, double h, int & zone)
{
Vec3d p1p (globp1, locpoint);
auto& gi = geominfo.GetPGI(1);
auto n = geo.GetNormal(gi.trignum, locpoint, &gi);
auto p1p = locpoint - p1;
plainpoint(0) = (p1p * ex) / h;
plainpoint(1) = (p1p * ey) / h;
// p1p = locpoint - globp1;
p1p /= h;
plainpoint[0] = p1p * ex;
plainpoint[1] = p1p * ey;
if(n*ez < 0)
zone = -1;
else
zone = 0;
}
@ -155,12 +174,9 @@ namespace netgen
PointGeomInfo & gi,
double h)
{
Vec3d p1p;
gi.trignum = 1;
p1p = plainpoint[0] * ex + plainpoint[1] * ey;
p1p *= h;
locpoint = globp1 + p1p;
locpoint = p1 + (h*plainpoint(0)) * ex + (h* plainpoint(1)) * ey;
if (!geo.ProjectPointGI(gi.trignum, locpoint, gi))
gi = geo.ProjectPoint(gi.trignum, locpoint);
return 0;
}
@ -853,6 +869,7 @@ namespace netgen
for (int i = oldnp+1; i <= plainpoints.Size(); i++)
{
Point<3> locp;
upgeominfo.Elem(i) = *blgeominfo1;
int err =
TransformFromPlain (plainpoints.Elem(i), locp,
upgeominfo.Elem(i), h);

View File

@ -41,12 +41,16 @@ class Meshing2
///
double maxarea;
Vec3d ex, ey;
Point3d globp1;
Vec3d ex, ey, ez;
Point<3> p1, p2;
const NetgenGeometry& geo;
public:
///
DLL_HEADER Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox);
DLL_HEADER Meshing2 (const NetgenGeometry& geo,
const MeshingParameters & mp,
const Box<3> & aboundingbox);
///
DLL_HEADER virtual ~Meshing2 ();

View File

@ -2801,7 +2801,9 @@ namespace netgen
<< " elementorder = " << elementorder << endl
<< " quad = " << quad << endl
<< " inverttets = " << inverttets << endl
<< " inverttrigs = " << inverttrigs << endl;
<< " inverttrigs = " << inverttrigs << endl
<< "closeedge enabled = " << closeedgefac.has_value() << endl
<< "closeedgefac = " << closeedgefac.value_or(0.) << endl;
}
/*

View File

@ -271,6 +271,9 @@ namespace netgen
void DoArchive (Archive & ar) { ar & i; }
};
inline void SetInvalid (SurfaceElementIndex & id) { id = -1; }
inline bool IsInvalid (SurfaceElementIndex & id) { return id == -1; }
inline istream & operator>> (istream & ist, SurfaceElementIndex & pi)
{
int i; ist >> i; pi = i; return ist;
@ -1231,7 +1234,7 @@ namespace netgen
// P .. plot, pause
// c .. combine
**/
string optimize2d = "smsmsmSmSmSm";
string optimize2d = "smcmSmcmSmcm";
/// number of 2d optimization steps
int optsteps2d = 3;
/// power of error (to approximate max err optimization)
@ -1256,6 +1259,8 @@ namespace netgen
double minh = 0.0;
/// file for meshsize
string meshsizefilename = "";
/// restrict h based on close edges
optional<double> closeedgefac = nullopt;
/// start surfacemeshing from everywhere in surface
bool startinsurface = false;
/// check overlapping surfaces (debug)

View File

@ -3,7 +3,7 @@
namespace netgen
{
int printmessage_importance = 5;
int printmessage_importance = 3;
int printwarnings = 1;
int printerrors = 1;
int printdots = 1;

View File

@ -321,35 +321,35 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
;
py::class_<Element2d>(m, "Element2D")
.def(py::init ([](int index, py::list vertices)
.def(py::init ([](int index, std::vector<PointIndex> vertices)
{
Element2d * newel = nullptr;
if (py::len(vertices) == 3)
if (vertices.size() == 3)
{
newel = new Element2d(TRIG);
for (int i = 0; i < 3; i++)
(*newel)[i] = py::extract<PointIndex>(vertices[i])();
(*newel)[i] = vertices[i];
newel->SetIndex(index);
}
else if (py::len(vertices) == 4)
else if (vertices.size() == 4)
{
newel = new Element2d(QUAD);
for (int i = 0; i < 4; i++)
(*newel)[i] = py::extract<PointIndex>(vertices[i])();
(*newel)[i] = vertices[i];
newel->SetIndex(index);
}
else if (py::len(vertices) == 6)
else if (vertices.size() == 6)
{
newel = new Element2d(TRIG6);
for(int i = 0; i<6; i++)
(*newel)[i] = py::extract<PointIndex>(vertices[i])();
(*newel)[i] = vertices[i];
newel->SetIndex(index);
}
else if (py::len(vertices) == 8)
else if (vertices.size() == 8)
{
newel = new Element2d(QUAD8);
for(int i = 0; i<8; i++)
(*newel)[i] = py::extract<PointIndex>(vertices[i])();
(*newel)[i] = vertices[i];
newel->SetIndex(index);
}
else
@ -801,6 +801,20 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
return self.Compress ();
} ,py::call_guard<py::gil_scoped_release>())
.def ("AddRegion", [] (Mesh & self, string name, int dim) -> int
{
auto & regionnames = self.GetRegionNamesCD(self.GetDimension()-dim);
regionnames.Append (new string(name));
int idx = regionnames.Size();
if (dim == 2)
{
FaceDescriptor fd;
fd.SetBCName(regionnames.Last());
fd.SetBCProperty(idx);
self.AddFaceDescriptor(fd);
}
return idx;
}, py::arg("name"), py::arg("dim"))
.def ("SetBCName", &Mesh::SetBCName)
.def ("GetBCName", FunctionPointer([](Mesh & self, int bc)->string
@ -983,6 +997,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
res["tet"] = py::make_tuple( values[2], values[3] );
return res;
}, py::arg("badelement_limit")=175.0)
.def ("CalcTotalBadness", &Mesh::CalcTotalBad)
.def ("GetQualityHistogram", &Mesh::GetQualityHistogram)
;
m.def("ImportMesh", [](const string& filename)

View File

@ -66,7 +66,7 @@ optimize3d: str = "cmdmustm"
optsteps3d: int = 3
Number of 3d optimization steps.
optimize2d: str = "smsmsmSmSmSm"
optimize2d: str = "smcmSmcmSmcm"
2d optimization strategy:
s .. swap, opt 6 lines/node
S .. swap, optimal elements
@ -176,6 +176,8 @@ inline void CreateMPfromKwargs(MeshingParameters& mp, py::kwargs kwargs, bool th
mp.parallel_meshing = py::cast<bool>(kwargs.attr("pop")("parallel_meshing"));
if(kwargs.contains("nthreads"))
mp.nthreads = py::cast<int>(kwargs.attr("pop")("nthreads"));
if(kwargs.contains("closeedgefac"))
mp.closeedgefac = py::cast<optional<double>>(kwargs.attr("pop")("closeedgefac"));
if(kwargs.size())
{
if(throw_if_not_all_parsed)

View File

@ -147,7 +147,7 @@ namespace netgen
{
pointset[pinew] = true;
Point<3> pnew;
PointBetween (mesh.Point (el[0]),
geo.PointBetweenEdge(mesh.Point (el[0]),
mesh.Point (el[1]), 0.5,
el.surfnr1, el.surfnr2,
el.epgeominfo[0], el.epgeominfo[1],
@ -216,7 +216,7 @@ namespace netgen
Point<3> pb;
PointGeomInfo pgi;
PointBetween (mesh.Point (pi1),
geo.PointBetween(mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]),
@ -307,7 +307,7 @@ namespace netgen
else
{
Point<3> pb;
PointBetween (mesh.Point (pi1),
geo.PointBetween(mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]),

View File

@ -100,7 +100,7 @@ namespace netgen
{
Point<3> pb;
EdgePointGeomInfo ngi;
PointBetween (mesh.Point (el[0]),
geo.PointBetweenEdge(mesh.Point (el[0]),
mesh.Point (el[1]), 0.5,
el.surfnr1, el.surfnr2,
el.epgeominfo[0], el.epgeominfo[1],
@ -184,7 +184,7 @@ namespace netgen
{
Point<3> pb;
PointGeomInfo newgi;
PointBetween (mesh.Point (pi1),
geo.PointBetween(mesh.Point (pi1),
mesh.Point (pi2), 0.5,
mesh.GetFaceDescriptor(el.GetIndex ()).SurfNr(),
el.GeomInfoPi (betw[j][0]+1),

View File

@ -15,14 +15,14 @@ namespace netgen
if(surfaceindex[i] >= 0)
{
*dest[i] = *from[i];
ProjectPoint(surfaceindex[i],*dest[i]);
geo.ProjectPoint(surfaceindex[i],*dest[i]);
}
}
}
void MeshOptimize2d :: ImproveVolumeMesh (Mesh & mesh)
void MeshOptimize2d :: ImproveVolumeMesh ()
{
if (!faceindex)
@ -31,7 +31,7 @@ namespace netgen
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
{
ImproveVolumeMesh (mesh);
ImproveVolumeMesh ();
if (multithread.terminate)
throw NgException ("Meshing stopped");
}
@ -229,7 +229,7 @@ namespace netgen
//cout << "origp " << origp << " newp " << mesh[pi];
ngi = gi1;
moveisok = (ProjectPointGI (surfi, mesh[pi], ngi) != 0);
moveisok = (geo.ProjectPointGI(surfi, mesh[pi], ngi) != 0);
//cout << " projected " << mesh[pi] << endl;

View File

@ -205,22 +205,20 @@ namespace netgen
class Opti2SurfaceMinFunction : public MinFunction
{
const Mesh & mesh;
Opti2dLocalData & ld;
const NetgenGeometry& geo;
public:
Opti2SurfaceMinFunction (const Mesh & amesh,
Opti2dLocalData & ald)
: mesh(amesh), ld(ald)
: ld(ald), geo(*amesh.GetGeometry())
{ } ;
virtual double Func (const Vector & x) const
{
Vec<3> n;
double badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
Point<3> pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
for (int j = 0; j < ld.locelements.Size(); j++)
@ -355,13 +353,13 @@ namespace netgen
// static int timer = NgProfiler::CreateTimer ("opti2surface - funcgrad");
// NgProfiler::RegionTimer reg (timer);
Vec<3> n, vgrad;
Vec<3> vgrad;
Point<3> pp1;
vgrad = 0;
double badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
// meshthis -> ProjectPoint (surfi, pp1);
@ -410,13 +408,13 @@ namespace netgen
// static int timer = NgProfiler::CreateTimer ("opti2surface - funcderiv");
// NgProfiler::RegionTimer reg (timer);
Vec<3> n, vgrad;
Vec<3> vgrad;
Point<3> pp1;
vgrad = 0;
double badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
for (int j = 0; j < ld.locelements.Size(); j++)
@ -474,11 +472,12 @@ namespace netgen
{
const Mesh & mesh;
Opti2dLocalData & ld;
const NetgenGeometry& geo;
public:
Opti2EdgeMinFunction (const Mesh & amesh,
Opti2dLocalData & ald)
: mesh(amesh), ld(ald) { } ;
: mesh(amesh), ld(ald), geo(*amesh.GetGeometry()) { } ;
virtual double FuncGrad (const Vector & x, Vector & g) const;
virtual double Func (const Vector & x) const;
@ -493,7 +492,7 @@ namespace netgen
double Opti2EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
{
int j, rot;
Vec<3> n1, n2, v1, v2, e1, e2, vgrad;
Vec<3> v1, v2, e1, e2, vgrad;
Point<3> pp1;
Vec<2> g1;
double badness, hbadness;
@ -502,7 +501,7 @@ namespace netgen
badness = 0;
pp1 = ld.sp1 + x(0) * ld.t1;
ld.meshthis -> ProjectPoint2 (ld.surfi, ld.surfi2, pp1);
geo.ProjectPointEdge(ld.surfi, ld.surfi2, pp1);
for (j = 0; j < ld.locelements.Size(); j++)
{
@ -526,8 +525,8 @@ namespace netgen
vgrad += g1(0) * e1 + g1(1) * e2;
}
ld.meshthis -> GetNormalVector (ld.surfi, pp1, n1);
ld.meshthis -> GetNormalVector (ld.surfi2, pp1, n2);
auto n1 = geo.GetNormal(ld.surfi, pp1);
auto n2 = geo.GetNormal(ld.surfi2, pp1);
v1 = Cross (n1, n2);
v1.Normalize();
@ -544,11 +543,12 @@ namespace netgen
{
const Mesh & mesh;
Opti2dLocalData & ld;
const NetgenGeometry& geo;
public:
Opti2SurfaceMinFunctionJacobian (const Mesh & amesh,
Opti2dLocalData & ald)
: mesh(amesh), ld(ald)
: mesh(amesh), ld(ald), geo(*amesh.GetGeometry())
{ } ;
virtual double FuncGrad (const Vector & x, Vector & g) const;
virtual double FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const;
@ -569,7 +569,7 @@ namespace netgen
// from 2d:
int lpi, gpi;
Vec<3> n, vgrad;
Vec<3> vgrad;
Point<3> pp1;
Vec<2> g1, vdir;
double badness, hbad, hderiv;
@ -577,7 +577,7 @@ namespace netgen
vgrad = 0;
badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
// auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
@ -641,7 +641,7 @@ namespace netgen
// from 2d:
int j, k, lpi, gpi;
Vec<3> n, vgrad;
Vec<3> vgrad;
Point<3> pp1;
Vec<2> g1, vdir;
double badness, hbad, hderiv;
@ -649,8 +649,6 @@ namespace netgen
vgrad = 0;
badness = 0;
ld.meshthis -> GetNormalVector (ld.surfi, ld.sp1, ld.gi1, n);
// pp1 = sp1;
// pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
@ -690,60 +688,25 @@ namespace netgen
return badness;
}
MeshOptimize2d dummy;
MeshOptimize2d :: MeshOptimize2d ()
void MeshOptimize2d :: ImproveMesh (const MeshingParameters & mp)
{
SetFaceIndex (0);
SetImproveEdges (0);
SetMetricWeight (0);
SetWriteStatus (1);
}
static Timer timer("MeshSmoothing 2D"); RegionTimer reg (timer);
void MeshOptimize2d :: SelectSurfaceOfPoint (const Point<3> & p,
const PointGeomInfo & gi)
{
;
}
void MeshOptimize2d :: ImproveMesh (Mesh & mesh, const MeshingParameters & mp)
{
if (!faceindex)
{
PrintMessage (3, "Smoothing");
for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
{
ImproveMesh (mesh, mp);
if (multithread.terminate)
throw NgException ("Meshing stopped");
}
faceindex = 0;
return;
}
static Timer timer("MeshSmoothing 2D");
// static int timer1 = NgProfiler::CreateTimer ("MeshSmoothing 2D start");
// static int timer2 = NgProfiler::CreateTimer ("MeshSmoothing 2D - BFGS");
RegionTimer reg (timer);
// NgProfiler::StartTimer (timer1);
CheckMeshApproximation (mesh);
Opti2dLocalData ld;
int ncolors;
Array<int> colors;
bool mixed = false;
auto elementsonpoint = mesh.CreatePoint2SurfaceElementTable( faceindex );
NgArray<MeshPoint, PointIndex::BASE> savepoints(mesh.GetNP());
Table<PointIndex> color_table;
if(faceindex)
{
Array<SurfaceElementIndex> seia;
mesh.GetSurfaceElementsOfFace (faceindex, seia);
bool mixed = 0;
for (auto sei : seia)
if (mesh[sei].GetNP() != 3)
{
@ -751,13 +714,7 @@ namespace netgen
break;
}
Vector x(2);
NgArray<MeshPoint, PointIndex::BASE> savepoints(mesh.GetNP());
ld.uselocalh = mp.uselocalh;
NgArray<int, PointIndex::BASE> compress(mesh.GetNP());
Array<int, PointIndex> compress(mesh.GetNP());
NgArray<PointIndex> icompress;
for (int i = 0; i < seia.Size(); i++)
{
@ -775,58 +732,52 @@ namespace netgen
icompress.Append(el[j]);
}
}
NgArray<int> cnta(icompress.Size());
cnta = 0;
for (int i = 0; i < seia.Size(); i++)
const auto & getDofs = [&] (int i)
{
const Element2d & el = mesh[seia[i]];
for (int j = 0; j < el.GetNP(); j++)
cnta[compress[el[j]]]++;
}
TABLE<SurfaceElementIndex> elementsonpoint(cnta);
for (int i = 0; i < seia.Size(); i++)
return elementsonpoint[icompress[i]];
};
colors.SetSize(icompress.Size());
ncolors = ngcore::ComputeColoring( colors, mesh.GetNSE(), getDofs );
TableCreator<PointIndex> creator(ncolors);
for ( ; !creator.Done(); creator++)
ParallelForRange( Range(colors), [&](auto myrange)
{
const Element2d & el = mesh[seia[i]];
for (int j = 0; j < el.GetNP(); j++)
elementsonpoint.Add (compress[el[j]], seia[i]);
for(auto i : myrange)
creator.Add(colors[i], icompress[i]);
});
color_table = creator.MoveTable();
}
/*
NgArray<int, PointIndex::BASE> nelementsonpoint(mesh.GetNP());
nelementsonpoint = 0;
for (int i = 0; i < seia.Size(); i++)
else
{
const Element2d & el = mesh[seia[i]];
for (int j = 0; j < el.GetNP(); j++)
nelementsonpoint[el[j]]++;
}
TABLE<SurfaceElementIndex,PointIndex::BASE> elementsonpoint(nelementsonpoint);
for (int i = 0; i < seia.Size(); i++)
for (auto & se : mesh.SurfaceElements())
if (se.GetNP() != 3)
{
const Element2d & el = mesh[seia[i]];
for (int j = 0; j < el.GetNP(); j++)
elementsonpoint.Add (el[j], seia[i]);
mixed = true;
break;
}
*/
const auto & getDofs = [&] (int i)
{
return elementsonpoint[i+PointIndex::BASE];
};
colors.SetSize(mesh.GetNP());
ncolors = ngcore::ComputeColoring( colors, mesh.GetNSE(), getDofs );
TableCreator<PointIndex> creator(ncolors);
for ( ; !creator.Done(); creator++)
ParallelForRange( Range(colors), [&](auto myrange)
{
for(auto i : myrange)
creator.Add(colors[i], PointIndex(i+PointIndex::BASE));
});
ld.loch = mp.maxh;
ld.locmetricweight = metricweight;
ld.meshthis = this;
Opti2SurfaceMinFunction surfminf(mesh, ld);
Opti2EdgeMinFunction edgeminf(mesh, ld);
Opti2SurfaceMinFunctionJacobian surfminfj(mesh, ld);
OptiParameters par;
par.maxit_linsearch = 8;
par.maxit_bfgs = 5;
color_table = creator.MoveTable();
}
/*
int i, j, k;
@ -897,27 +848,6 @@ namespace netgen
*/
bool printeddot = 0;
char plotchar = '.';
int modplot = 1;
if (mesh.GetNP() > 1000)
{
plotchar = '+';
modplot = 100;
}
if (mesh.GetNP() > 10000)
{
plotchar = 'o';
modplot = 1000;
}
if (mesh.GetNP() > 100000)
{
plotchar = 'O';
modplot = 10000;
}
int cnt = 0;
// NgProfiler::StopTimer (timer1);
/*
@ -927,28 +857,39 @@ namespace netgen
static Timer tloop("MeshSmooting 2D - loop");
tloop.Start();
for (int hi = 0; hi < icompress.Size(); hi++)
for (auto icolor : Range(color_table))
{
PointIndex pi = icompress[hi];
if (multithread.terminate)
break;
ParallelForRange( Range(color_table[icolor].Size()), [&](auto myrange)
{
Opti2dLocalData ld;
ld.uselocalh = mp.uselocalh;
ld.loch = mp.maxh;
ld.locmetricweight = metricweight;
ld.meshthis = this;
Opti2SurfaceMinFunction surfminf(mesh, ld);
Opti2SurfaceMinFunctionJacobian surfminfj(mesh, ld);
MinFunction & minfunc = mixed ? static_cast<MinFunction&>(surfminfj) : surfminf;
OptiParameters par;
par.maxit_linsearch = 8;
par.maxit_bfgs = 5;
for (auto i : myrange)
{
PointIndex pi = color_table[icolor][i];
if (mesh[pi].Type() == SURFACEPOINT)
{
if (multithread.terminate)
throw NgException ("Meshing stopped");
return;
cnt++;
if (cnt % modplot == 0 && writestatus)
{
printeddot = 1;
PrintDot (plotchar);
}
// if (elementsonpoint[pi].Size() == 0) continue;
if (elementsonpoint[hi].Size() == 0) continue;
if (elementsonpoint[pi].Size() == 0) continue;
ld.sp1 = mesh[pi];
// Element2d & hel = mesh[elementsonpoint[pi][0]];
Element2d & hel = mesh[elementsonpoint[hi][0]];
Element2d & hel = mesh[elementsonpoint[pi][0]];
int hpi = 0;
for (int j = 1; j <= hel.GetNP(); j++)
@ -959,7 +900,7 @@ namespace netgen
}
ld.gi1 = hel.GeomInfoPi(hpi);
SelectSurfaceOfPoint (ld.sp1, ld.gi1);
// SelectSurfaceOfPoint (ld.sp1, ld.gi1);
ld.locelements.SetSize(0);
ld.locrots.SetSize (0);
@ -967,9 +908,9 @@ namespace netgen
ld.loc_pnts2.SetSize (0);
ld.loc_pnts3.SetSize (0);
for (int j = 0; j < elementsonpoint[hi].Size(); j++)
for (int j = 0; j < elementsonpoint[pi].Size(); j++)
{
SurfaceElementIndex sei = elementsonpoint[hi][j];
SurfaceElementIndex sei = elementsonpoint[pi][j];
const Element2d & bel = mesh[sei];
ld.surfi = mesh.GetFaceDescriptor(bel.GetIndex()).SurfNr();
@ -992,11 +933,13 @@ namespace netgen
}
GetNormalVector (ld.surfi, ld.sp1, ld.gi1, ld.normal);
ld.normal = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1);
ld.t1 = ld.normal.GetNormal ();
ld.t2 = Cross (ld.normal, ld.t1);
// save points, and project to tangential plane
if(mixed)
{
// save points, and project to tangential plane (only for optimization with Opti2SurfaceMinFunctionJacobian in mixed element meshes)
for (int j = 0; j < ld.locelements.Size(); j++)
{
const Element2d & el = mesh[ld.locelements[j]];
@ -1014,28 +957,25 @@ namespace netgen
mesh[hhpi] -= lam * ld.normal;
}
}
}
Vector x(2);
x = 0;
par.typx = 0.3*ld.lochs[0];
// NgProfiler::StartTimer (timer2);
if (mixed)
{
BFGS (x, surfminfj, par, 1e-6);
}
else
{
BFGS (x, surfminf, par, 1e-6);
}
BFGS (x, minfunc, par, 1e-6);
// NgProfiler::StopTimer (timer2);
Point3d origp = mesh[pi];
auto origp = mesh[pi];
int loci = 1;
double fact = 1;
int moveisok = 0;
if(mixed)
{
// restore other points
for (int j = 0; j < ld.locelements.Size(); j++)
{
@ -1046,6 +986,7 @@ namespace netgen
if (hhpi != pi) mesh[hhpi] = savepoints[hhpi];
}
}
}
//optimizer loop (if whole distance is not possible, move only a bit!!!!)
@ -1070,7 +1011,7 @@ namespace netgen
PointGeomInfo ngi;
ngi = ld.gi1;
moveisok = ProjectPointGI (ld.surfi, mesh[pi], ngi);
moveisok = geo.ProjectPointGI(ld.surfi, mesh[pi], ngi);
// point lies on same chart in stlsurface
if (moveisok)
@ -1080,28 +1021,17 @@ namespace netgen
}
else
{
mesh[pi] = Point<3> (origp);
mesh[pi] = origp;
}
}
}
}
}, mixed ? 1 : ngcore::TasksPerThread(4)); // mixed element smoothing not parallel yet
}
tloop.Stop();
if (printeddot)
PrintDot ('\n');
CheckMeshApproximation (mesh);
mesh.SetNextTimeStamp();
}
void MeshOptimize2d :: GetNormalVector(INDEX /* surfind */, const Point<3> & p, Vec<3> & nv) const
{
nv = Vec<3> (0, 0, 1);
}
void MeshOptimize2d :: GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & gi, Vec<3> & n) const
{
GetNormalVector (surfind, p, n);
}
}

View File

@ -930,46 +930,6 @@ double Opti3EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
double CalcTotalBad (const Mesh::T_POINTS & points,
const Array<Element> & elements,
const MeshingParameters & mp)
{
static Timer t("CalcTotalBad"); RegionTimer reg(t);
static constexpr int n_classes = 20;
double sum = 0;
tets_in_qualclass.SetSize(n_classes);
tets_in_qualclass = 0;
ParallelForRange( IntRange(elements.Size()), [&] (auto myrange) {
double local_sum = 0.0;
double teterrpow = mp.opterrpow;
std::array<int,n_classes> classes_local{};
for (auto i : myrange)
{
double elbad = pow (max2(CalcBad (points, elements[i], 0, mp),1e-10),
1/teterrpow);
int qualclass = int (n_classes / elbad + 1);
if (qualclass < 1) qualclass = 1;
if (qualclass > n_classes) qualclass = n_classes;
classes_local[qualclass-1]++;
local_sum += elbad;
}
AtomicAdd(sum, local_sum);
for (auto i : Range(n_classes))
AsAtomic(tets_in_qualclass[i]) += classes_local[i];
});
return sum;
}
int WrongOrientation (const Mesh::T_POINTS & points, const Element & el)
{
const Point3d & p1 = points[el.PNum(1)];
@ -1204,7 +1164,7 @@ void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal)
}
const char * savetask = multithread.task;
multithread.task = "Smooth Mesh";
multithread.task = "Optimize Volume: Smooth Mesh";
TABLE<INDEX> surfelementsonpoint(points.Size());
@ -1383,7 +1343,7 @@ void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL g
if (goal == OPT_QUALITY)
{
double bad1 = CalcTotalBad (points, volelements, mp);
double bad1 = CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
PrintMessage (5, "Total badness = ", bad1);
}
@ -1438,7 +1398,7 @@ void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL g
const char * savetask = multithread.task;
multithread.task = "Smooth Mesh";
multithread.task = "Optimize Volume: Smooth Mesh";
for (PointIndex pi : points.Range())
if ( (*this)[pi].Type() == INNERPOINT )
@ -1485,7 +1445,7 @@ void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL g
if (goal == OPT_QUALITY)
{
double bad1 = CalcTotalBad (points, volelements, mp);
double bad1 = CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
PrintMessage (5, "Total badness = ", bad1);
}
@ -1536,7 +1496,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
if (goal == OPT_QUALITY)
{
double bad1 = CalcTotalBad (points, volelements, mp);
double bad1 = CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
PrintMessage (5, "Total badness = ", bad1);
}
@ -1564,7 +1524,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
}
const char * savetask = multithread.task;
multithread.task = "Smooth Mesh";
multithread.task = "Optimize Volume: Smooth Mesh";
topt.Start();
int counter = 0;
@ -1631,7 +1591,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
if (goal == OPT_QUALITY)
{
double bad1 = CalcTotalBad (points, volelements, mp);
double bad1 = CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
PrintMessage (5, "Total badness = ", bad1);
}
@ -1699,7 +1659,7 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
const char * savetask = multithread.task;
multithread.task = "Smooth Mesh Jacobian";
multithread.task = "Optimize Volume: Smooth Mesh Jacobian";
// for (PointIndex pi = points.Begin(); i < points.End(); pi++)
for (PointIndex pi : points.Range())
@ -1855,7 +1815,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
const char * savetask = multithread.task;
multithread.task = "Smooth Mesh Jacobian";
multithread.task = "Optimize Volume: Smooth Mesh Jacobian";
// for (PointIndex pi = points.Begin(); pi <= points.End(); pi++)
for (PointIndex pi : points.Range())

View File

@ -276,8 +276,8 @@ namespace netgen
double oldlamedge,oldlamface;
MeshOptimize2d * optimizer2d = refinement.Get2dOptimizer();
if(!optimizer2d)
auto geo = mesh.GetGeometry();
if(!geo)
{
cerr << "No 2D Optimizer!" << endl;
return;
@ -382,8 +382,15 @@ namespace netgen
for (int i = 1; i <= np; i++)
*can.Elem(i) = mesh.Point(i);
if(optimizer2d)
optimizer2d->ProjectBoundaryPoints(surfaceindex,can,should);
if(geo)
for(int i=0; i<surfaceindex.Size(); i++)
{
if(surfaceindex[i] >= 0)
{
*should[i] = *can[i];
geo->ProjectPoint(surfaceindex[i],*should[i]);
}
}
}

View File

@ -12,6 +12,11 @@ target_link_libraries(occ PUBLIC ngcore)
if(NOT WIN32)
target_link_libraries( occ PRIVATE ${OCC_LIBRARIES} ${PYTHON_LIBRARIES})
if(USE_OCC AND APPLE)
# Link AppKit in case OCE was built as static libraries
find_library(AppKit AppKit)
target_link_libraries( occ PRIVATE ${AppKit} )
endif(USE_OCC AND APPLE)
install( TARGETS occ ${NG_INSTALL_DIR})
if (USE_GUI)
target_link_libraries( occvis PUBLIC occ )

View File

@ -206,7 +206,7 @@ static void PutInBounds (const TopoDS_Face& F,
Handle (Geom_Surface) S = BRep_Tool::Surface(F,L);
if (S->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
S = (*(Handle_Geom_RectangularTrimmedSurface*)&S)->BasisSurface();
S = Handle(Geom_RectangularTrimmedSurface)::DownCast(S)->BasisSurface();
}
if (!S->IsUPeriodic() && !S->IsVPeriodic())
return;
@ -702,7 +702,7 @@ TopTools_MapOfShape& Partition_Inter3d::TouchedFaces()
//purpose :
//=======================================================================
Handle_BRepAlgo_AsDes Partition_Inter3d::AsDes() const
Handle(BRepAlgo_AsDes) Partition_Inter3d::AsDes() const
{
return myAsDes;
}
@ -829,7 +829,7 @@ TopoDS_Vertex Partition_Inter3d::ReplaceSameDomainV(const TopoDS_Vertex& V,
//purpose :
//=======================================================================
Handle_BRepAlgo_AsDes Partition_Inter3d::SectionEdgesAD() const
Handle(BRepAlgo_AsDes) Partition_Inter3d::SectionEdgesAD() const
{
return mySectionEdgesAD;
}

View File

@ -96,13 +96,13 @@ public:
void FacesPartition(const TopoDS_Face& F1,const TopoDS_Face& F2) ;
Standard_Boolean IsDone(const TopoDS_Face& F1,const TopoDS_Face& F2) const;
TopTools_MapOfShape& TouchedFaces() ;
Handle_BRepAlgo_AsDes AsDes() const;
Handle(BRepAlgo_AsDes) AsDes() const;
TopTools_MapOfShape& NewEdges() ;
Standard_Boolean HasSameDomainF(const TopoDS_Shape& F) const;
Standard_Boolean IsSameDomainF(const TopoDS_Shape& F1,const TopoDS_Shape& F2) const;
const TopTools_ListOfShape& SameDomain(const TopoDS_Face& F) const;
TopoDS_Vertex ReplaceSameDomainV(const TopoDS_Vertex& V,const TopoDS_Edge& E) const;
Handle_BRepAlgo_AsDes SectionEdgesAD() const;
Handle(BRepAlgo_AsDes) SectionEdgesAD() const;
Standard_Boolean IsSectionEdge(const TopoDS_Edge& E) const;
Standard_Boolean HasSectionEdge(const TopoDS_Face& F) const;
Standard_Boolean IsSplitOn(const TopoDS_Edge& NewE,const TopoDS_Edge& OldE,const TopoDS_Face& F) const;
@ -134,11 +134,11 @@ private:
// Fields PRIVATE
//
Handle_BRepAlgo_AsDes myAsDes;
Handle(BRepAlgo_AsDes) myAsDes;
TopTools_DataMapOfShapeListOfShape myDone;
TopTools_MapOfShape myTouched;
TopTools_MapOfShape myNewEdges;
Handle_BRepAlgo_AsDes mySectionEdgesAD;
Handle(BRepAlgo_AsDes) mySectionEdgesAD;
TopTools_DataMapOfShapeListOfShape mySameDomainFM;
TopTools_DataMapOfShapeShape mySameDomainVM;

View File

@ -143,7 +143,7 @@ private:
TopTools_DataMapOfShapeShape myFaceShapeMap;
TopTools_DataMapOfShapeShape myInternalFaces;
TopTools_DataMapOfShapeShape myIntNotClFaces;
Handle_BRepAlgo_AsDes myAsDes;
Handle(BRepAlgo_AsDes) myAsDes;
BRepAlgo_Image myImagesFaces;
BRepAlgo_Image myImagesEdges;
BRepAlgo_Image myImageShape;

View File

@ -69,7 +69,7 @@ namespace netgen
void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h,
BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h,
const MeshingParameters & mparam)
{
int ls = -1;
@ -112,41 +112,41 @@ namespace netgen
{
double curvature = 0;
prop->SetParameters (parmid.X(), parmid.Y());
if (!prop->IsCurvatureDefined())
prop2->SetParameters (parmid.X(), parmid.Y());
if (!prop2->IsCurvatureDefined())
{
(*testout) << "curvature not defined!" << endl;
return;
}
curvature = max(fabs(prop->MinCurvature()),
fabs(prop->MaxCurvature()));
curvature = max(fabs(prop2->MinCurvature()),
fabs(prop2->MaxCurvature()));
prop->SetParameters (par0.X(), par0.Y());
if (!prop->IsCurvatureDefined())
prop2->SetParameters (par0.X(), par0.Y());
if (!prop2->IsCurvatureDefined())
{
(*testout) << "curvature not defined!" << endl;
return;
}
curvature = max(curvature,max(fabs(prop->MinCurvature()),
fabs(prop->MaxCurvature())));
curvature = max(curvature,max(fabs(prop2->MinCurvature()),
fabs(prop2->MaxCurvature())));
prop->SetParameters (par1.X(), par1.Y());
if (!prop->IsCurvatureDefined())
prop2->SetParameters (par1.X(), par1.Y());
if (!prop2->IsCurvatureDefined())
{
(*testout) << "curvature not defined!" << endl;
return;
}
curvature = max(curvature,max(fabs(prop->MinCurvature()),
fabs(prop->MaxCurvature())));
curvature = max(curvature,max(fabs(prop2->MinCurvature()),
fabs(prop2->MaxCurvature())));
prop->SetParameters (par2.X(), par2.Y());
if (!prop->IsCurvatureDefined())
prop2->SetParameters (par2.X(), par2.Y());
if (!prop2->IsCurvatureDefined())
{
(*testout) << "curvature not defined!" << endl;
return;
}
curvature = max(curvature,max(fabs(prop->MinCurvature()),
fabs(prop->MaxCurvature())));
curvature = max(curvature,max(fabs(prop2->MinCurvature()),
fabs(prop2->MaxCurvature())));
//(*testout) << "curvature " << curvature << endl;
@ -165,7 +165,7 @@ namespace netgen
return;
if (h > 30) return;
// if (h > 30) return;
}
if (h < maxside && depth < 10)
@ -181,20 +181,20 @@ namespace netgen
if(ls == 0)
{
pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y()));
RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam);
}
else if(ls == 1)
{
pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y()));
RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam);
}
else if(ls == 2)
{
pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y()));
RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam);
}
}
@ -306,7 +306,7 @@ namespace netgen
void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam)
void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam)
{
static Timer t("OCCFindEdges"); RegionTimer r(t);
static Timer tsearch("OCCFindEdges - search point");
@ -601,8 +601,8 @@ namespace netgen
void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh,
MeshingParameters & mparam)
void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh,
const MeshingParameters & mparam)
{
static Timer t("OCCMeshSurface"); RegionTimer r(t);
@ -647,11 +647,11 @@ namespace netgen
Box<3> bb = geom.GetBoundingBox();
// int projecttype = PLANESPACE;
int projecttype = PLANESPACE;
static Timer tinit("init");
tinit.Start();
Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam);
Meshing2OCCSurfaces meshing(geom, TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam);
tinit.Stop();
@ -796,7 +796,6 @@ namespace netgen
// Philippose - 15/01/2009
double maxh = geom.face_maxh[k-1];
//double maxh = mparam.maxh;
mparam.checkoverlap = 0;
// int noldpoints = mesh->GetNP();
int noldsurfel = mesh.GetNSE();
@ -809,9 +808,13 @@ namespace netgen
MESHING2_RESULT res;
// TODO: check overlap not correctly working here
MeshingParameters mparam_without_overlap = mparam;
mparam_without_overlap.checkoverlap = false;
try {
static Timer t("GenerateMesh"); RegionTimer reg(t);
res = meshing.GenerateMesh (mesh, mparam, maxh, k);
res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k);
}
catch (SingularMatrixException)
@ -916,7 +919,7 @@ namespace netgen
}
void OCCOptimizeSurface(OCCGeometry & geom, Mesh & mesh,
MeshingParameters & mparam)
const MeshingParameters & mparam)
{
const char * savetask = multithread.task;
multithread.task = "Optimizing surface";
@ -941,41 +944,41 @@ namespace netgen
if (multithread.terminate) return;
{
MeshOptimize2dOCCSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2));
meshopt.EdgeSwapping (i > mparam.optsteps2d/2);
}
if (multithread.terminate) return;
{
MeshOptimize2dOCCSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mesh, mparam);
meshopt.ImproveMesh (mparam);
}
{
MeshOptimize2dOCCSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.CombineImprove (mesh);
meshopt.CombineImprove ();
}
if (multithread.terminate) return;
{
MeshOptimize2dOCCSurfaces meshopt(geom);
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mesh, mparam);
meshopt.ImproveMesh (mparam);
}
}
@ -991,9 +994,11 @@ namespace netgen
void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh,
void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh,
const MeshingParameters & mparam, const OCCParameters& occparam)
{
static Timer t1("OCCSetLocalMeshSize");
RegionTimer regt(t1);
mesh.SetGlobalH (mparam.maxh);
mesh.SetMinimalH (mparam.minh);
@ -1029,11 +1034,15 @@ namespace netgen
multithread.task = "Setting local mesh size (elements per edge)";
// setting elements per edge
// Philippose - 23/01/2009
// Find all the parent faces of a given edge
// and limit the mesh size of the edge based on the
// mesh size limit of the face
TopTools_IndexedDataMapOfShapeListOfShape edge_face_map;
edge_face_map.Clear();
TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map);
// setting elements per edge
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
TopoDS_Edge e = TopoDS::Edge (geom.emap(i));
@ -1137,7 +1146,9 @@ namespace netgen
}
BRepAdaptor_Surface sf(face, Standard_True);
BRepLProp_SLProps prop(sf, 2, 1e-5);
// one prop for evaluating and one for derivatives
BRepLProp_SLProps prop(sf, 0, 1e-5);
BRepLProp_SLProps prop2(sf, 2, 1e-5);
int ntriangles = triangulation -> NbTriangles();
for (int j = 1; j <= ntriangles; j++)
@ -1158,14 +1169,14 @@ namespace netgen
//maxside = max (maxside, p[1].Distance(p[2]));
//cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;
RestrictHTriangle (par[0], par[1], par[2], &prop, mesh, 0, 0, mparam);
RestrictHTriangle (par[0], par[1], par[2], &prop, &prop2, mesh, 0, 0, mparam);
//cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
}
}
// setting close edges
if (occparam.resthcloseedgeenable)
if (mparam.closeedgefac.has_value())
{
multithread.task = "Setting local mesh size (close edges)";
@ -1250,7 +1261,7 @@ namespace netgen
mindist = min (mindist, line.Dist(lines[num]));
}
mindist /= (occparam.resthcloseedgefac + VSMALL);
mindist /= (*mparam.closeedgefac + VSMALL);
if (mindist < 1e-3 * bb.Diam())
{
@ -1273,197 +1284,6 @@ namespace netgen
mesh.LoadLocalMeshSize (mparam.meshsizefilename);
}
int OCCGenerateMesh (OCCGeometry & geom, shared_ptr<Mesh> & mesh, MeshingParameters & mparam,
const OCCParameters& occparam)
{
multithread.percent = 0;
if (mparam.perfstepsstart <= MESHCONST_ANALYSE)
{
if(mesh.get() == nullptr)
mesh = make_shared<Mesh>();
mesh->geomtype = Mesh::GEOM_OCC;
OCCSetLocalMeshSize(geom,*mesh, mparam, occparam);
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_ANALYSE)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_MESHEDGES)
{
OCCFindEdges (geom, *mesh, mparam);
/*
cout << "Removing redundant points" << endl;
int i, j;
int np = mesh->GetNP();
NgArray<int> equalto;
equalto.SetSize (np);
equalto = 0;
for (i = 1; i <= np; i++)
{
for (j = i+1; j <= np; j++)
{
if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12))
equalto[j-1] = i;
}
}
for (i = 1; i <= np; i++)
if (equalto[i-1])
{
cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl;
for (j = 1; j <= mesh->GetNSeg(); j++)
{
Segment & seg = mesh->LineSegment(j);
if (seg[0] == i) seg[0] = equalto[i-1];
if (seg[1] == i) seg[1] = equalto[i-1];
}
}
cout << "Removing degenerated segments" << endl;
for (j = 1; j <= mesh->GetNSeg(); j++)
{
Segment & seg = mesh->LineSegment(j);
if (seg[0] == seg[1])
{
mesh->DeleteSegment(j);
cout << "Deleting Segment " << j << endl;
}
}
mesh->Compress();
*/
/*
for (int i = 1; i <= geom.fmap.Extent(); i++)
{
Handle(Geom_Surface) hf1 =
BRep_Tool::Surface(TopoDS::Face(geom.fmap(i)));
for (int j = i+1; j <= geom.fmap.Extent(); j++)
{
Handle(Geom_Surface) hf2 =
BRep_Tool::Surface(TopoDS::Face(geom.fmap(j)));
if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl;
}
}
*/
#ifdef LOG_STREAM
(*logout) << "Edges meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
{
OCCMeshSurface (geom, *mesh, mparam);
if (multithread.terminate) return TCL_OK;
#ifdef LOG_STREAM
(*logout) << "Surfaces meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
#ifdef STAT_STREAM
(*statout) << mesh->GetNSeg() << " & "
<< mesh->GetNSE() << " & - &"
<< GetTime() << " & " << endl;
#endif
// MeshQuality2d (*mesh);
mesh->CalcSurfacesOfNode();
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHSURFACE)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_OPTSURFACE)
{
OCCOptimizeSurface(geom, *mesh, mparam);
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
{
multithread.task = "Volume meshing";
MESHING3_RESULT res = MeshVolume (mparam, *mesh);
if (res != MESHING3_OK) return TCL_ERROR;
if (multithread.terminate) return TCL_OK;
RemoveIllegalElements (*mesh);
if (multithread.terminate) return TCL_OK;
MeshQuality3d (*mesh);
#ifdef STAT_STREAM
(*statout) << GetTime() << " & ";
#endif
#ifdef LOG_STREAM
(*logout) << "Volume meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHVOLUME)
return TCL_OK;
if (mparam.perfstepsstart <= MESHCONST_OPTVOLUME)
{
multithread.task = "Volume optimization";
OptimizeVolume (mparam, *mesh);
if (multithread.terminate) return TCL_OK;
#ifdef STAT_STREAM
(*statout) << GetTime() << " & "
<< mesh->GetNE() << " & "
<< mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl;
#endif
#ifdef LOG_STREAM
(*logout) << "Volume optimized" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
// cout << "Optimization complete" << endl;
}
/*
(*testout) << "NP: " << mesh->GetNP() << endl;
for (int i = 1; i <= mesh->GetNP(); i++)
(*testout) << mesh->Point(i) << endl;
(*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl;
for (int i = 1; i <= mesh->GetNSeg(); i++)
(*testout) << mesh->LineSegment(i) << endl;
*/
for (int i = 0; i < mesh->GetNDomains(); i++)
if (geom.snames.Size())
mesh->SetMaterial (i+1, geom.snames[i]);
return TCL_OK;
}
}
#endif

View File

@ -74,6 +74,30 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
}
}
void OCCGeometry :: Analyse(Mesh& mesh,
const MeshingParameters& mparam) const
{
OCCSetLocalMeshSize(*this, mesh, mparam, occparam);
}
void OCCGeometry :: FindEdges(Mesh& mesh,
const MeshingParameters& mparam) const
{
OCCFindEdges(*this, mesh, mparam);
}
void OCCGeometry :: MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) const
{
OCCMeshSurface(*this, mesh, mparam);
}
void OCCGeometry :: FinalizeMesh(Mesh& mesh) const
{
for (int i = 0; i < mesh.GetNDomains(); i++)
if (snames.Size())
mesh.SetMaterial (i+1, snames[i]);
}
void OCCGeometry :: PrintNrShapes ()
{
@ -1010,10 +1034,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
SetCenter();
}
void OCCGeometry :: Project (int surfi, Point<3> & p) const
PointGeomInfo OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const
{
static int cnt = 0;
if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
@ -1027,13 +1048,55 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
suval.Coord( u, v);
pnt = thesurf->Value( u, v );
PointGeomInfo gi;
gi.trignum = surfi;
gi.u = u;
gi.v = v;
p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
return gi;
}
bool OCCGeometry :: ProjectPointGI(int surfind, Point<3>& p, PointGeomInfo& gi) const
{
double u = gi.u;
double v = gi.v;
Point<3> hp = p;
if (FastProject (surfind, hp, u, v))
{
p = hp;
return 1;
}
ProjectPoint (surfind, p);
return CalcPointGeomInfo (surfind, gi, p);
}
void OCCGeometry :: ProjectPointEdge(int surfind, INDEX surfind2,
Point<3> & p, EdgePointGeomInfo* gi) const
{
TopExp_Explorer exp0, exp1;
bool done = false;
Handle(Geom_Curve) c;
for (exp0.Init(fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next())
for (exp1.Init(fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next())
{
if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current())))
{
done = true;
double s0, s1;
c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1);
}
}
gp_Pnt pnt(p(0), p(1), p(2));
GeomAPI_ProjectPointOnCurve proj(pnt, c);
pnt = proj.NearestPoint();
p(0) = pnt.X();
p(1) = pnt.Y();
p(2) = pnt.Z();
}
bool OCCGeometry :: FastProject (int surfi, Point<3> & ap, double& u, double& v) const
{
@ -1091,7 +1154,147 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
return true;
}
Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* geominfo) const
{
if(geominfo)
{
gp_Pnt pnt;
gp_Vec du, dv;
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
occface->D1(geominfo->u,geominfo->v,pnt,du,dv);
auto n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
Vec<3>(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
return n;
}
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) );
suval.Coord( u, v);
pnt = occface->Value( u, v );
gp_Vec du, dv;
occface->D1(u,v,pnt,du,dv);
/*
if (!occface->IsCNu (1) || !occface->IsCNv (1))
(*testout) << "SurfOpt: Differentiation FAIL" << endl;
*/
auto n = Cross (Vec3d(du.X(), du.Y(), du.Z()),
Vec3d(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
return n;
}
bool OCCGeometry :: CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const
{
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return 0;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) );
suval.Coord( u, v);
//pnt = occface->Value( u, v );
gi.u = u;
gi.v = v;
return true;
}
void OCCGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi > 0)
{
double u = gi1.u+secpoint*(gi2.u-gi1.u);
double v = gi1.v+secpoint*(gi2.v-gi1.v);
auto savept = hnewp;
if (!FastProject(surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2))
{
// cout << "Fast projection to surface fails! Using OCC projection" << endl;
hnewp = savept;
ProjectPoint(surfi, hnewp);
}
newgi.trignum = 1;
newgi.u = u;
newgi.v = v;
}
newp = hnewp;
}
void OCCGeometry :: PointBetweenEdge(const Point<3> & p1,
const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
double s0, s1;
Point<3> hnewp = p1+secpoint*(p2-p1);
gp_Pnt pnt(hnewp(0), hnewp(1), hnewp(2));
GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(emap(ap1.edgenr)), s0, s1));
pnt = proj.NearestPoint();
hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
newp = hnewp;
newgi = ap1;
};
// void OCCGeometry :: WriteOCC_STL(char * filename)
@ -1109,6 +1312,10 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
void LoadOCCInto(OCCGeometry* occgeo, const char* filename)
{
static Timer timer_all("LoadOCC"); RegionTimer rtall(timer_all);
static Timer timer_readfile("LoadOCC-ReadFile");
static Timer timer_transfer("LoadOCC-Transfer");
static Timer timer_getnames("LoadOCC-get names");
// Initiate a dummy XCAF Application to handle the STEP XCAF Document
static Handle_XCAFApp_Application dummy_app = XCAFApp_Application::GetApplication();
@ -1125,19 +1332,23 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
}
dummy_app->NewDocument ("STEP-XCAF",step_doc);
timer_readfile.Start();
STEPCAFControl_Reader reader;
// Enable transfer of colours
reader.SetColorMode(Standard_True);
reader.SetNameMode(Standard_True);
Standard_Integer stat = reader.ReadFile((char*)filename);
timer_readfile.Stop();
timer_transfer.Start();
if(stat != IFSelect_RetDone)
{
throw NgException("Couldn't load OCC geometry");
}
reader.Transfer(step_doc);
timer_transfer.Stop();
// Read in the shape(s) and the colours present in the STEP File
Handle_XCAFDoc_ShapeTool step_shape_contents = XCAFDoc_DocumentTool::ShapeTool(step_doc->Main());
@ -1175,6 +1386,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
occgeo->snames.Append(name);
TopExp_Explorer exp0,exp1;
timer_getnames.Start();
for (exp0.Init(occgeo->shape, TopAbs_FACE); exp0.More(); exp0.Next())
{
TopoDS_Face face = TopoDS::Face(exp0.Current());
@ -1182,13 +1394,14 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
if (name == string(""))
snprintf(name, 50, "bc_%zu", occgeo->fnames.Size());
occgeo->fnames.Append(name);
for (exp1.Init(face, TopAbs_EDGE); exp1.More(); exp1.Next())
{
TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
STEP_GetEntityName(edge,&reader,name);
occgeo->enames.Append(name);
}
// for (exp1.Init(face, TopAbs_EDGE); exp1.More(); exp1.Next())
// {
// TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
// STEP_GetEntityName(edge,&reader,name);
// occgeo->enames.Append(name);
// }
}
timer_getnames.Stop();
// Gerhard BEGIN
// cout << "Solid Names: "<<endl;
// for (int i=0;i<occgeo->snames.Size();i++)
@ -1681,21 +1894,9 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
return false;
}
const Refinement & OCCGeometry :: GetRefinement () const
{
return * new OCCRefinementSurfaces (*this);
}
void OCCParameters :: Print(ostream & ost) const
{
ost << "OCC Parameters:" << endl
<< "close edges: " << resthcloseedgeenable
<< ", fac = " << resthcloseedgefac << endl
<< "minimum edge length: " << resthminedgelenenable
<< ", min len = " << resthminedgelen << endl;
}
@ -1703,10 +1904,10 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a
DLL_HEADER extern OCCParameters occparam;
OCCParameters occparam;
int OCCGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
{
return OCCGenerateMesh (*this, mesh, mparam, occparam);
}
// int OCCGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
// {
// return OCCGenerateMesh (*this, mesh, mparam, occparam);
// }
}

View File

@ -73,7 +73,6 @@
#include "ShapeUpgrade_ShellSewing.hxx"
#include "ShapeFix_Shape.hxx"
#include "ShapeFix_Wireframe.hxx"
#include "BRepMesh.hxx"
#include "BRepMesh_IncrementalMesh.hxx"
#include "BRepBndLib.hxx"
#include "Bnd_Box.hxx"
@ -184,19 +183,40 @@ namespace netgen
return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00;
}
class DLL_HEADER OCCParameters
{
public:
/// Factor for meshing close edges, moved to meshingparameters
// double resthcloseedgefac = 2.;
/// Enable / Disable detection of close edges
// int resthcloseedgeenable = true;
/// Minimum edge length to be used for dividing edges to mesh points
double resthminedgelen = 0.001;
/// Enable / Disable use of the minimum edge length (by default use 1e-4)
int resthminedgelenenable = true;
/*!
Dump all the OpenCascade specific meshing parameters
to console
*/
void Print (ostream & ost) const;
};
class OCCGeometry : public NetgenGeometry
{
Point<3> center;
OCCParameters occparam;
public:
TopoDS_Shape shape;
TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap;
NgArray<bool> fsingular, esingular, vsingular;
Box<3> boundingbox;
NgArray<string> fnames, enames, snames;
NgArray<string> fnames, /*enames,*/ snames;
// Philippose - 29/01/2009
// OpenCascade XDE Support
// XCAF Handle to make the face colours available to the rest of
@ -204,7 +224,7 @@ namespace netgen
Handle_XCAFDoc_ColorTool face_colours;
mutable int changed;
NgArray<int> facemeshstatus;
mutable NgArray<int> facemeshstatus;
// Philippose - 15/01/2009
// Maximum mesh size for a given face
@ -241,10 +261,42 @@ namespace netgen
vmap.Clear();
}
Mesh::GEOM_TYPE GetGeomType() const override
{ return Mesh::GEOM_OCC; }
DLL_HEADER virtual void Save (string filename) const;
void SetOCCParameters(const OCCParameters& par)
{ occparam = par; }
void DoArchive(Archive& ar);
void Analyse(Mesh& mesh,
const MeshingParameters& mparam) const override;
void FindEdges(Mesh& mesh,
const MeshingParameters& mparam) const override;
void MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) const override;
void FinalizeMesh(Mesh& mesh) const override;
DLL_HEADER void Save (string filename) const override;
void DoArchive(Archive& ar) override;
PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override;
void ProjectPointEdge (int surfind, int surfind2, Point<3> & p,
EdgePointGeomInfo* gi = nullptr) const override;
bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override;
bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
DLL_HEADER void BuildFMap();
@ -265,9 +317,6 @@ namespace netgen
Point<3> Center() const
{ return center; }
void Project (int surfi, Point<3> & p) const;
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
OCCSurface GetSurface (int surfi)
{
cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;
@ -392,37 +441,9 @@ namespace netgen
// void WriteOCC_STL(char * filename);
DLL_HEADER virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
DLL_HEADER virtual const Refinement & GetRefinement () const;
};
class DLL_HEADER OCCParameters
{
public:
/// Factor for meshing close edges
double resthcloseedgefac = 2.;
/// Enable / Disable detection of close edges
int resthcloseedgeenable = true;
/// Minimum edge length to be used for dividing edges to mesh points
double resthminedgelen = 0.001;
/// Enable / Disable use of the minimum edge length (by default use 1e-4)
int resthminedgelenenable = true;
/*!
Dump all the OpenCascade specific meshing parameters
to console
*/
void Print (ostream & ost) const;
// DLL_HEADER virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
private:
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
};
@ -435,17 +456,14 @@ namespace netgen
// Philippose - 31.09.2009
// External access to the mesh generation functions within the OCC
// subsystem (Not sure if this is the best way to implement this....!!)
DLL_HEADER extern int OCCGenerateMesh (OCCGeometry & occgeometry, shared_ptr<Mesh> & mesh,
MeshingParameters & mparam, const OCCParameters& occparam);
DLL_HEADER extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam,
DLL_HEADER extern void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam,
const OCCParameters& occparam);
DLL_HEADER extern void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, MeshingParameters & mparam);
DLL_HEADER extern void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
DLL_HEADER extern void OCCOptimizeSurface (OCCGeometry & geom, Mesh & mesh, MeshingParameters & mparam);
DLL_HEADER extern void OCCOptimizeSurface (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
DLL_HEADER extern void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
DLL_HEADER extern void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
}
#endif

View File

@ -343,6 +343,8 @@ namespace netgen
PointGeomInfo & gi,
double h)
{
static Timer t("FromPlane"); RegionTimer reg(t);
if (projecttype == PLANESPACE)
{
// cout << "2d : " << pplane << endl;
@ -366,12 +368,76 @@ namespace netgen
void OCCSurface :: Project (Point<3> & p, PointGeomInfo & gi)
void OCCSurface :: Project (Point<3> & ap, PointGeomInfo & gi)
{
static Timer t("OccSurface::Project"); RegionTimer reg(t);
static Timer t2("OccSurface::Project actural");
// try Newton's method ...
gp_Pnt p(ap(0), ap(1), ap(2));
double u = gi.u;
double v = gi.v;
gp_Pnt x = occface->Value (u,v);
if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return;
gp_Vec du, dv;
occface->D1(u,v,x,du,dv);
int count = 0;
gp_Pnt xold;
gp_Vec n;
double det, lambda, mu;
do
{
count++;
n = du^dv;
det = Det3 (n.X(), du.X(), dv.X(),
n.Y(), du.Y(), dv.Y(),
n.Z(), du.Z(), dv.Z());
if (det < 1e-15)
break;
lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
n.Y(), p.Y()-x.Y(), dv.Y(),
n.Z(), p.Z()-x.Z(), dv.Z())/det;
mu = Det3 (n.X(), du.X(), p.X()-x.X(),
n.Y(), du.Y(), p.Y()-x.Y(),
n.Z(), du.Z(), p.Z()-x.Z())/det;
u += lambda;
v += mu;
xold = x;
occface->D1(u,v,x,du,dv);
if (xold.SquareDistance(x) < sqr(PROJECTION_TOLERANCE))
{
ap = Point<3> (x.X(), x.Y(), x.Z());
gi.u = u;
gi.v = v;
return;
}
}
while (count < 20);
// Newton did not converge, use OCC projection
// static int cnt = 0;
// if (cnt++ % 1000 == 0) cout << "********************************************** OCCSurfce :: Project, cnt = " << cnt << endl;
gp_Pnt pnt(p(0), p(1), p(2));
gp_Pnt pnt = p; // (p(0), p(1), p(2));
//(*testout) << "pnt = " << pnt.X() << ", " << pnt.Y() << ", " << pnt.Z() << endl;
@ -406,27 +472,32 @@ namespace netgen
proj.LowerDistanceParameters (gi.u, gi.v);
*/
double u,v;
// double u,v;
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( topods_face ) );
auto toltool = BRep_Tool::Tolerance( topods_face );
// gp_Pnt2d suval = su->ValueOfUV ( pnt, toltool);
t2.Start();
gp_Pnt2d suval = su->NextValueOfUV (gp_Pnt2d(u,v), pnt, toltool);
t2.Stop();
suval.Coord( u, v);
pnt = occface->Value( u, v );
//(*testout) << "pnt(proj) = " << pnt.X() << ", " << pnt.Y() << ", " << pnt.Z() << endl;
gi.u = u;
gi.v = v;
gi.trignum = 1;
p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
ap = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
}
Meshing2OCCSurfaces :: Meshing2OCCSurfaces (const TopoDS_Shape & asurf,
Meshing2OCCSurfaces :: Meshing2OCCSurfaces (const NetgenGeometry& geo,
const TopoDS_Shape & asurf,
const Box<3> & abb, int aprojecttype,
const MeshingParameters & mparam)
: Meshing2(mparam, Box<3>(abb.PMin(), abb.PMax())), surface(TopoDS::Face(asurf), aprojecttype)
: Meshing2(geo, mparam, Box<3>(abb.PMin(), abb.PMax())),
surface(TopoDS::Face(asurf), aprojecttype)
{
;
}
@ -463,188 +534,6 @@ namespace netgen
return gh;
}
MeshOptimize2dOCCSurfaces :: MeshOptimize2dOCCSurfaces (const OCCGeometry & ageometry)
: MeshOptimize2d(), geometry(ageometry)
{
;
}
void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
geometry.Project (surfind, p);
}
int MeshOptimize2dOCCSurfaces :: ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
{
double u = gi.u;
double v = gi.v;
Point<3> hp = p;
if (geometry.FastProject (surfind, hp, u, v))
{
p = hp;
return 1;
}
ProjectPoint (surfind, p);
return CalcPointGeomInfo (surfind, gi, p);
}
void MeshOptimize2dOCCSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2,
Point<3> & p) const
{
TopExp_Explorer exp0, exp1;
bool done = false;
Handle(Geom_Curve) c;
for (exp0.Init(geometry.fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next())
for (exp1.Init(geometry.fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next())
{
if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current())))
{
done = true;
double s0, s1;
c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1);
}
}
gp_Pnt pnt(p(0), p(1), p(2));
GeomAPI_ProjectPointOnCurve proj(pnt, c);
pnt = proj.NearestPoint();
p(0) = pnt.X();
p(1) = pnt.Y();
p(2) = pnt.Z();
}
void MeshOptimize2dOCCSurfaces ::
GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & geominfo, Vec<3> & n) const
{
gp_Pnt pnt;
gp_Vec du, dv;
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
occface->D1(geominfo.u,geominfo.v,pnt,du,dv);
n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
Vec<3>(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (geometry.fmap(surfind).Orientation() == TopAbs_REVERSED) n = -1*n;
// GetNormalVector (surfind, p, n);
}
void MeshOptimize2dOCCSurfaces ::
GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const
{
// static int cnt = 0;
// if (cnt++ % 1000 == 0) cout << "GetNV cnt = " << cnt << endl;
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(geometry.fmap(surfind)) ) );
suval.Coord( u, v);
pnt = occface->Value( u, v );
gp_Vec du, dv;
occface->D1(u,v,pnt,du,dv);
/*
if (!occface->IsCNu (1) || !occface->IsCNv (1))
(*testout) << "SurfOpt: Differentiation FAIL" << endl;
*/
n = Cross (Vec3d(du.X(), du.Y(), du.Z()),
Vec3d(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (geometry.fmap(surfind).Orientation() == TopAbs_REVERSED) n = -1*n;
}
int MeshOptimize2dOCCSurfaces ::
CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const
{
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return 0;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(geometry.fmap(surfind)) ) );
suval.Coord( u, v);
//pnt = occface->Value( u, v );
gi.u = u;
gi.v = v;
return 1;
}
OCCRefinementSurfaces :: OCCRefinementSurfaces (const OCCGeometry & ageometry)
: Refinement(), geometry(ageometry)
{
;
}
OCCRefinementSurfaces :: ~OCCRefinementSurfaces ()
{
;
}
/*
inline double Det3 (double a00, double a01, double a02,
double a10, double a11, double a12,
@ -703,76 +592,6 @@ namespace netgen
return true;
}
*/
void OCCRefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi > 0)
{
double u = gi1.u+secpoint*(gi2.u-gi1.u);
double v = gi1.v+secpoint*(gi2.v-gi1.v);
auto savept = hnewp;
if (!geometry.FastProject (surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2))
{
// cout << "Fast projection to surface fails! Using OCC projection" << endl;
hnewp = savept;
geometry.Project (surfi, hnewp);
}
newgi.trignum = 1;
newgi.u = u;
newgi.v = v;
}
newp = hnewp;
}
void OCCRefinementSurfaces ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
double s0, s1;
Point<3> hnewp = p1+secpoint*(p2-p1);
gp_Pnt pnt(hnewp(0), hnewp(1), hnewp(2));
GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(geometry.emap(ap1.edgenr)), s0, s1));
pnt = proj.NearestPoint();
hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
newp = hnewp;
newgi = ap1;
};
void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi) const
{
if (surfi > 0)
geometry.Project (surfi, p);
};
void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi) const
{
if (surfi > 0)
if (!geometry.FastProject (surfi, p, gi.u, gi.v))
{
cout << "Fast projection to surface fails! Using OCC projection" << endl;
geometry.Project (surfi, p);
}
};
}

View File

@ -113,7 +113,8 @@ class Meshing2OCCSurfaces : public Meshing2
public:
///
Meshing2OCCSurfaces (const TopoDS_Shape & asurf, const Box<3> & aboundingbox,
Meshing2OCCSurfaces (const NetgenGeometry& geo,
const TopoDS_Shape & asurf, const Box<3> & aboundingbox,
int aprojecttype, const MeshingParameters & mparam);
///
@ -141,64 +142,8 @@ protected:
};
///
class MeshOptimize2dOCCSurfaces : public MeshOptimize2d
{
///
const OCCGeometry & geometry;
public:
///
MeshOptimize2dOCCSurfaces (const OCCGeometry & ageometry);
///
virtual void ProjectPoint (INDEX surfind, Point<3> & p) const;
///
virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const;
///
virtual int ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const;
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const;
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & gi, Vec<3> & n) const;
virtual int CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const;
};
class OCCGeometry;
class DLL_HEADER OCCRefinementSurfaces : public Refinement
{
const OCCGeometry & geometry;
public:
OCCRefinementSurfaces (const OCCGeometry & ageometry);
virtual ~OCCRefinementSurfaces ();
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi) const override;
};
#endif

View File

@ -45,14 +45,12 @@ namespace netgen
virtual void SetParameters (Tcl_Interp * interp)
{
occparam.resthcloseedgefac =
atof (Tcl_GetVar (interp, "::stloptions.resthcloseedgefac", 0));
occparam.resthcloseedgeenable =
atoi (Tcl_GetVar (interp, "::stloptions.resthcloseedgeenable", 0));
occparam.resthminedgelen =
atof (Tcl_GetVar (interp, "::stloptions.resthminedgelen", 0));
occparam.resthminedgelenenable =
atoi (Tcl_GetVar (interp, "::stloptions.resthminedgelenenable", 0));
if(auto geo = dynamic_pointer_cast<OCCGeometry>(ng_geometry); geo)
geo->SetOCCParameters(occparam);
}
};

View File

@ -7,6 +7,7 @@
#include <meshing.hpp>
#include <occgeom.hpp>
#include <Standard_Version.hxx>
using namespace netgen;
@ -30,17 +31,6 @@ minedgelen: Optional[float] = 0.001
void CreateOCCParametersFromKwargs(OCCParameters& occparam, py::dict kwargs)
{
if(kwargs.contains("closeedgefac"))
{
auto val = kwargs.attr("pop")("closeedgefac");
if(val.is_none())
occparam.resthcloseedgeenable = false;
else
{
occparam.resthcloseedgefac = py::cast<double>(val);
occparam.resthcloseedgeenable = true;
}
}
if(kwargs.contains("minedgelen"))
{
auto val = kwargs.attr("pop")("minedgelen");
@ -57,6 +47,7 @@ void CreateOCCParametersFromKwargs(OCCParameters& occparam, py::dict kwargs)
DLL_HEADER void ExportNgOCC(py::module &m)
{
m.attr("occ_version") = OCC_VERSION_COMPLETE;
py::class_<OCCGeometry, shared_ptr<OCCGeometry>, NetgenGeometry> (m, "OCCGeometry", R"raw_string(Use LoadOCCGeometry to load the geometry from a *.step file.)raw_string")
.def(py::init<>())
.def(py::init([] (const string& filename)
@ -183,11 +174,12 @@ DLL_HEADER void ExportNgOCC(py::module &m)
CreateOCCParametersFromKwargs(occparam, kwargs);
CreateMPfromKwargs(mp, kwargs);
}
geo->SetOCCParameters(occparam);
auto mesh = make_shared<Mesh>();
SetGlobalMesh(mesh);
mesh->SetGeometry(geo);
geo->GenerateMesh(mesh, mp);
SetGlobalMesh(mesh);
ng_geometry = geo;
OCCGenerateMesh(*geo, mesh, mp, occparam);
return mesh;
}, py::arg("mp") = nullptr,
py::call_guard<py::gil_scoped_release>(),

View File

@ -22,7 +22,7 @@ static void STLFindEdges (STLGeometry & geom, Mesh & mesh,
// mark edge points:
//int ngp = geom.GetNP();
geom.RestrictLocalH(mesh, h, stlparam);
geom.RestrictLocalH(mesh, h, stlparam, mparam);
PushStatusF("Mesh Lines");
@ -299,15 +299,14 @@ int STLSurfaceMeshing (STLGeometry & geom, class Mesh & mesh, const MeshingParam
geom.SetMarkedTrig(seg.geominfo[1].trignum,1);
}
MeshOptimizeSTLSurface optmesh(geom);
MeshOptimize2d optmesh(mesh);
optmesh.SetFaceIndex (0);
optmesh.SetImproveEdges (0);
optmesh.SetMetricWeight (0);
mesh.CalcSurfacesOfNode();
optmesh.EdgeSwapping (mesh, 0);
mesh.CalcSurfacesOfNode();
optmesh.ImproveMesh (mesh, mparam);
optmesh.EdgeSwapping (0);
optmesh.ImproveMesh (mparam);
}
mesh.Compress();
@ -827,7 +826,7 @@ void STLSurfaceOptimization (STLGeometry & geom,
{
PrintFnStart("optimize STL Surface");
MeshOptimizeSTLSurface optmesh(geom);
MeshOptimize2d optmesh(mesh);
optmesh.SetFaceIndex (0);
optmesh.SetImproveEdges (0);
@ -848,25 +847,41 @@ void STLSurfaceOptimization (STLGeometry & geom,
{
case 's':
{
optmesh.EdgeSwapping (mesh, 0);
optmesh.EdgeSwapping(0);
break;
}
case 'S':
{
optmesh.EdgeSwapping (mesh, 1);
optmesh.EdgeSwapping(1);
break;
}
case 'm':
{
optmesh.ImproveMesh(mesh, mparam);
optmesh.ImproveMesh(mparam);
break;
}
case 'c':
{
optmesh.CombineImprove (mesh);
optmesh.CombineImprove();
break;
}
}
// while(mesh.CheckOverlappingBoundary())
// {
// for(const auto & el : mesh.SurfaceElements())
// {
// if(el.BadElement())
// {
// cout << "Restrict localh at el nr " << el << endl;
// for(const auto& p : el.PNums())
// {
// const auto& pnt = mesh[p];
// mesh.RestrictLocalH(pnt, 0.5*mesh.GetH(pnt));
// }
// }
// }
// optmesh.SplitImprove();
// }
//(*testout) << "optimize, after, step = " << meshparam.optimize2d[j-1] << mesh.Point (3679) << endl;
}
@ -882,7 +897,7 @@ void STLSurfaceOptimization (STLGeometry & geom,
MeshingSTLSurface :: MeshingSTLSurface (STLGeometry & ageom,
const MeshingParameters & mp)
: Meshing2(mp, ageom.GetBoundingBox()), geom(ageom)
: Meshing2(ageom, mp, ageom.GetBoundingBox()), geom(ageom)
{
;
}
@ -1052,195 +1067,4 @@ double MeshingSTLSurface :: Area () const
return geom.Area();
}
MeshOptimizeSTLSurface :: MeshOptimizeSTLSurface (STLGeometry & ageom)
: MeshOptimize2d(), geom(ageom)
{
;
}
void MeshOptimizeSTLSurface :: SelectSurfaceOfPoint (const Point<3> & p,
const PointGeomInfo & gi)
{
// (*testout) << "sel char: " << gi.trignum << endl;
geom.SelectChartOfTriangle (gi.trignum);
// geom.SelectChartOfPoint (p);
}
void MeshOptimizeSTLSurface :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
if (!geom.Project (p))
{
PrintMessage(7,"project failed");
if (!geom.ProjectOnWholeSurface(p))
{
PrintMessage(7, "project on whole surface failed");
}
}
// geometry.GetSurface(surfind)->Project (p);
}
void MeshOptimizeSTLSurface :: ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const
{
/*
ProjectToEdge ( geometry.GetSurface(surfind),
geometry.GetSurface(surfind2), p);
*/
}
int MeshOptimizeSTLSurface :: CalcPointGeomInfo(PointGeomInfo& gi, const Point<3> & p3) const
{
Point<3> hp = p3;
gi.trignum = geom.Project (hp);
if (gi.trignum) return 1;
return 0;
}
void MeshOptimizeSTLSurface :: GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const
{
n = geom.GetChartNormalVector();
}
RefinementSTLGeometry :: RefinementSTLGeometry (const STLGeometry & ageom)
: Refinement(), geom(ageom)
{
;
}
RefinementSTLGeometry :: ~RefinementSTLGeometry ()
{
;
}
void RefinementSTLGeometry ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
/*
(*testout) << "surf-between: p1 = " << p1 << ", p2 = " << p2
<< ", gi = " << gi1 << " - " << gi2 << endl;
*/
if (gi1.trignum > 0)
{
// ((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
Point<3> np1 = newp;
Point<3> np2 = newp;
((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
int tn1 = geom.Project (np1);
((STLGeometry&)geom).SelectChartOfTriangle (gi2.trignum);
int tn2 = geom.Project (np2);
newgi.trignum = tn1; //urspruengliche version
newp = np1; //urspruengliche version
if (!newgi.trignum)
{ newgi.trignum = tn2; newp = np2; }
if (!newgi.trignum) newgi.trignum = gi1.trignum;
/*
if (tn1 != 0 && tn2 != 0 && ((STLGeometry&)geom).GetAngle(tn1,tn2) < M_PI*0.05) {
newgi.trignum = tn1;
newp = np1;
}
else
{
newp = ((STLGeometry&)geom).PointBetween(p1, gi1.trignum, p2, gi2.trignum);
tn1 = ((STLGeometry&)geom).Project(newp);
newgi.trignum = tn1;
if (!tn1)
{
newp = Center (p1, p2);
newgi.trignum = 0;
}
}
*/
}
else
{
// (*testout) << "WARNING: PointBetween got geominfo = 0" << endl;
newp = p1+secpoint*(p2-p1);
newgi.trignum = 0;
}
// (*testout) << "newp = " << newp << ", ngi = " << newgi << endl;
}
void RefinementSTLGeometry ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & gi1,
const EdgePointGeomInfo & gi2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
/*
(*testout) << "edge-between: p1 = " << p1 << ", p2 = " << p2
<< ", gi1,2 = " << gi1 << ", " << gi2 << endl;
*/
/*
newp = Center (p1, p2);
((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
newgi.trignum = geom.Project (newp);
*/
int hi;
newgi.dist = (1.0-secpoint) * gi1.dist + secpoint*gi2.dist;
newgi.edgenr = gi1.edgenr;
/*
(*testout) << "p1 = " << p1 << ", p2 = " << p2 << endl;
(*testout) << "refedge: " << gi1.edgenr
<< " d1 = " << gi1.dist << ", d2 = " << gi2.dist << endl;
*/
newp = geom.GetLine (gi1.edgenr)->GetPointInDist (geom.GetPoints(), newgi.dist, hi);
// (*testout) << "newp = " << newp << endl;
}
void RefinementSTLGeometry :: ProjectToSurface (Point<3> & p, int surfi) const
{
cout << "RefinementSTLGeometry :: ProjectToSurface not implemented!" << endl;
};
void RefinementSTLGeometry :: ProjectToSurface (Point<3> & p, int surfi,
PointGeomInfo & gi) const
{
((STLGeometry&)geom).SelectChartOfTriangle (gi.trignum);
gi.trignum = geom.Project (p);
// if (!gi.trignum)
// cout << "projectSTL failed" << endl;
};
}

View File

@ -63,59 +63,5 @@ protected:
double Area () const override;
};
///
class MeshOptimizeSTLSurface : public MeshOptimize2d
{
///
STLGeometry & geom;
public:
///
MeshOptimizeSTLSurface (STLGeometry & ageom);
///
virtual void SelectSurfaceOfPoint (const Point<3> & p,
const PointGeomInfo & gi);
///
virtual void ProjectPoint (INDEX surfind, Point<3> & p) const;
///
virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point<3> & p) const;
///
virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point<3> & p3) const;
///
virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const;
};
class RefinementSTLGeometry : public Refinement
{
const STLGeometry & geom;
public:
RefinementSTLGeometry (const STLGeometry & ageom);
virtual ~RefinementSTLGeometry ();
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi) const override;
virtual void ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi) const override;
};
#endif

View File

@ -52,7 +52,7 @@ void CreateSTLParametersFromKwargs(STLParameters& stlparam, py::dict kwargs)
stlparam.outerchartangle = py::cast<double>(kwargs.attr("pop")("outerchartangle"));
if(kwargs.contains("usesearchtree"))
stlparam.usesearchtree = py::cast<int>(kwargs.attr("pop")("usesearchtree"));
if(kwargs.contains("resthatlasfac"))
if(kwargs.contains("atlasfac"))
{
auto val = kwargs.attr("pop")("resthatlasfac");
if(val.is_none())
@ -65,9 +65,9 @@ void CreateSTLParametersFromKwargs(STLParameters& stlparam, py::dict kwargs)
}
if(kwargs.contains("atlasminh"))
stlparam.atlasminh = py::cast<double>(kwargs.attr("pop")("atlasminh"));
if(kwargs.contains("resthsurfcurvfac"))
if(kwargs.contains("surfcurvfac"))
{
auto val = kwargs.attr("pop")("resthsurfcurvfac");
auto val = kwargs.attr("pop")("surfcurvfac");
if(val.is_none())
stlparam.resthsurfcurvenable = false;
else
@ -76,9 +76,9 @@ void CreateSTLParametersFromKwargs(STLParameters& stlparam, py::dict kwargs)
stlparam.resthsurfcurvfac = py::cast<double>(val);
}
}
if(kwargs.contains("resthchartdistfac"))
if(kwargs.contains("chartdistfac"))
{
auto val = kwargs.attr("pop")("resthchartdistfac");
auto val = kwargs.attr("pop")("chartdistfac");
if(val.is_none())
stlparam.resthchartdistenable = false;
else
@ -87,20 +87,9 @@ void CreateSTLParametersFromKwargs(STLParameters& stlparam, py::dict kwargs)
stlparam.resthchartdistfac = py::cast<double>(val);
}
}
if(kwargs.contains("resthcloseedgefac"))
if(kwargs.contains("edgeanglefac"))
{
auto val = kwargs.attr("pop")("resthcloseedgefac");
if(val.is_none())
stlparam.resthcloseedgeenable = false;
else
{
stlparam.resthcloseedgeenable = true;
stlparam.resthcloseedgefac = py::cast<double>(val);
}
}
if(kwargs.contains("resthedgeanglefac"))
{
auto val = kwargs.attr("pop")("resthedgeanglefac");
auto val = kwargs.attr("pop")("edgeanglefac");
if(val.is_none())
stlparam.resthedgeangleenable = false;
else
@ -109,9 +98,9 @@ void CreateSTLParametersFromKwargs(STLParameters& stlparam, py::dict kwargs)
stlparam.resthedgeanglefac = py::cast<double>(val);
}
}
if(kwargs.contains("resthsurfmeshcurvfac"))
if(kwargs.contains("surfmeshcurvfac"))
{
auto val = kwargs.attr("pop")("resthsurfmeshcurvfac");
auto val = kwargs.attr("pop")("surfmeshcurvfac");
if(val.is_none())
stlparam.resthsurfmeshcurvenable = false;
else
@ -120,9 +109,9 @@ void CreateSTLParametersFromKwargs(STLParameters& stlparam, py::dict kwargs)
stlparam.resthsurfmeshcurvfac = py::cast<double>(val);
}
}
if(kwargs.contains("resthlinelengthfac"))
if(kwargs.contains("linelengthfac"))
{
auto val = kwargs.attr("pop")("resthlinelengthfac");
auto val = kwargs.attr("pop")("linelengthfac");
if(val.is_none())
stlparam.resthlinelengthenable = false;
else

View File

@ -44,7 +44,6 @@ void STLMeshing (STLGeometry & geom,
lineendpoints(), spiralpoints(), selectedmultiedge()
*/
{
ref = NULL;
edgedata = make_unique<STLEdgeDataList>(*this);
externaledges.SetSize(0);
Clear();
@ -66,7 +65,6 @@ STLGeometry :: ~STLGeometry()
{
// for (auto p : atlas) delete p;
// delete edgedata;
delete ref;
}
void STLGeometry :: Save (string filename) const
@ -102,17 +100,124 @@ int STLGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mp
return STLMeshingDummy (this, mesh, mparam, stlpar);
}
const Refinement & STLGeometry :: GetRefinement () const
Vec<3> STLGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const
{
delete ref;
ref = new RefinementSTLGeometry(*this);
// ref -> Set2dOptimizer(new MeshOptimizeSTLSurface(*this)); ??? copied from CSG
return *ref;
if(!gi)
throw Exception("STLGeometry::GetNormal without PointGeomInfo called");
return GetChart(GetChartNr(gi->trignum)).GetNormal();
}
bool STLGeometry :: CalcPointGeomInfo(int /*surfind*/, PointGeomInfo& gi, const Point<3> & p3) const
{
Point<3> hp = p3;
SelectChartOfTriangle(gi.trignum);
gi.trignum = Project (hp);
if (gi.trignum) return true;
return false;
}
bool STLGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
{
static std::mutex mutex_project_whole_surface;
int meshchart = GetChartNr(gi.trignum);
const STLChart& chart = GetChart(meshchart);
int trignum = chart.ProjectNormal(p);
if(trignum==0)
{
// non-thread-safe implementation
std::lock_guard<std::mutex> guard(mutex_project_whole_surface);
PrintMessage(7,"project failed");
SelectChartOfTriangle (gi.trignum); // needed because ProjectOnWholeSurface uses meshchartnv (the normal vector of selected chart)
trignum = ProjectOnWholeSurface(p);
if(trignum==0)
{
PrintMessage(7, "project on whole surface failed");
return false;
}
}
return true;
}
PointGeomInfo STLGeometry :: ProjectPoint (INDEX surfind, Point<3> & p) const
{
throw Exception("ProjectPoint without PointGeomInfo not implemented");
}
void STLGeometry ::
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
newp = p1+secpoint*(p2-p1);
/*
(*testout) << "surf-between: p1 = " << p1 << ", p2 = " << p2
<< ", gi = " << gi1 << " - " << gi2 << endl;
*/
if (gi1.trignum > 0)
{
// ((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
Point<3> np1 = newp;
Point<3> np2 = newp;
auto ngi1 = gi1;
auto ngi2 = gi2;
// SelectChartOfTriangle (gi1.trignum);
int tn1 = ProjectPointGI (surfi, np1, ngi1);
// SelectChartOfTriangle (gi2.trignum);
int tn2 = ProjectPointGI (surfi, np2, ngi2);
newgi.trignum = tn1; //urspruengliche version
newp = np1; //urspruengliche version
if (!newgi.trignum)
{ newgi.trignum = tn2; newp = np2; }
if (!newgi.trignum) newgi.trignum = gi1.trignum;
}
else
{
// (*testout) << "WARNING: PointBetween got geominfo = 0" << endl;
newp = p1+secpoint*(p2-p1);
newgi.trignum = 0;
}
}
void STLGeometry ::
PointBetweenEdge (const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & gi1,
const EdgePointGeomInfo & gi2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
/*
(*testout) << "edge-between: p1 = " << p1 << ", p2 = " << p2
<< ", gi1,2 = " << gi1 << ", " << gi2 << endl;
*/
/*
newp = Center (p1, p2);
((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
newgi.trignum = geom.Project (newp);
*/
int hi;
newgi.dist = (1.0-secpoint) * gi1.dist + secpoint*gi2.dist;
newgi.edgenr = gi1.edgenr;
/*
(*testout) << "p1 = " << p1 << ", p2 = " << p2 << endl;
(*testout) << "refedge: " << gi1.edgenr
<< " d1 = " << gi1.dist << ", d2 = " << gi2.dist << endl;
*/
newp = GetLine (gi1.edgenr)->GetPointInDist (GetPoints(), newgi.dist, hi);
// (*testout) << "newp = " << newp << endl;
}
void STLGeometry :: STLInfo(double* data)
{

View File

@ -148,7 +148,7 @@ namespace netgen
//for meshing and project:
NgArray<int> meshcharttrigs; //per trig: 1=belong to chart, 0 not
int meshchart;
mutable int meshchart;
NgArray<int> ha_points; // help array, np long, filled with 0
@ -159,12 +159,10 @@ namespace netgen
//transformation:
Vec<3> meshtrignv;
mutable Vec<3> meshtrignv;
Vec<3> ex, ey, ez;
Point<3> p1;
mutable class RefinementSTLGeometry * ref;
public:
int edgesfound;
int surfacemeshed;
@ -194,6 +192,23 @@ namespace netgen
virtual void Save (string filename) const override;
bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
PointGeomInfo ProjectPoint(INDEX surfind, Point<3> & p) const override;
bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const override;
void PointBetween(const Point<3> & p1, const Point<3> & p2,
double secpoint, int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
DLL_HEADER void STLInfo(double* data);
//stldoctor:
@ -419,7 +434,7 @@ namespace netgen
//
void DefineTangentialPlane(const Point<3> & ap1, const Point<3> & ap2, int trig);
//
void SelectChartOfTriangle (int trignum);
void SelectChartOfTriangle (int trignum) const;
//
void SelectChartOfPoint (const Point<3> & p);
//
@ -450,7 +465,7 @@ namespace netgen
int LineEndPointsSet() const {return lineendpoints.Size() == GetNP();}
void ClearLineEndPoints();
DLL_HEADER void RestrictLocalH(class Mesh & mesh, double gh, const STLParameters& stlparam);
DLL_HEADER void RestrictLocalH(class Mesh & mesh, double gh, const STLParameters& stlparam, const MeshingParameters& mparam);
void RestrictLocalHCurv(class Mesh & mesh, double gh, const STLParameters& stlparam);
void RestrictHChartDistOneChart(ChartId chartnum, NgArray<int>& acttrigs, class Mesh & mesh,
double gh, double fact, double minh, const STLParameters& stlparam);
@ -459,8 +474,6 @@ namespace netgen
int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam) override;
virtual const Refinement & GetRefinement () const override;
// Add additional Point to chart to close the surface and write the resulting stl to a file
DLL_HEADER void WriteChartToFile( ChartId chartnumber, string filename="chart.slb" );
};

View File

@ -392,7 +392,7 @@ void STLGeometry :: DefineTangentialPlane (const Point<3> & ap1, const Point<3>
}
void STLGeometry :: SelectChartOfTriangle (int trignum)
void STLGeometry :: SelectChartOfTriangle (int trignum) const
{
meshchart = GetChartNr(trignum);
meshtrignv = GetTriangle(trignum).Normal();
@ -814,7 +814,7 @@ void STLGeometry :: RestrictLocalHCurv(class Mesh & mesh, double gh, const STLPa
}
//restrict local h due to near edges and due to outer chart distance
void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh, const STLParameters& stlparam)
void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh, const STLParameters& stlparam, const MeshingParameters& mparam)
{
//bei jedem Dreieck alle Nachbardreiecke vergleichen, und, fallskein Kante dazwischen,
@ -921,12 +921,12 @@ void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh, const STLParame
PopStatus();
}
if (stlparam.resthcloseedgeenable)
if (mparam.closeedgefac.has_value())
{
PushStatusF("Restrict H due to close edges");
//geht nicht für spiralen!!!!!!!!!!!!!!!!!!
double disttohfact = sqr(10.0 / stlparam.resthcloseedgefac);
double disttohfact = sqr(10.0 / *mparam.closeedgefac);
int k,l;
double h1, h2, dist;
int rc = 0;

View File

@ -78,11 +78,6 @@ namespace netgen
stlparam.resthlinelengthenable =
atoi (Tcl_GetVar (interp, "::stloptions.resthlinelengthenable", 0));
stlparam.resthcloseedgefac =
atof (Tcl_GetVar (interp, "::stloptions.resthcloseedgefac", 0));
stlparam.resthcloseedgeenable =
atoi (Tcl_GetVar (interp, "::stloptions.resthcloseedgeenable", 0));
stlparam.resthedgeanglefac =
atof (Tcl_GetVar (interp, "::stloptions.resthedgeanglefac", 0));
stlparam.resthedgeangleenable =
@ -538,7 +533,7 @@ namespace netgen
mesh -> SetLocalH (stlgeometry->GetBoundingBox().PMin() - Vec3d(10, 10, 10),
stlgeometry->GetBoundingBox().PMax() + Vec3d(10, 10, 10),
mparam.grading);
stlgeometry -> RestrictLocalH(*mesh, mparam.maxh, stlparam);
stlgeometry -> RestrictLocalH(*mesh, mparam.maxh, stlparam, mparam);
if (stlparam.resthsurfmeshcurvenable)
mesh -> CalcLocalHFromSurfaceCurvature (mparam.grading,

View File

@ -1462,8 +1462,8 @@ STLParameters :: STLParameters()
resthchartdistenable = 1;
resthlinelengthfac = 0.5;
resthlinelengthenable = 1;
resthcloseedgefac = 1;
resthcloseedgeenable = 1;
// resthcloseedgefac = 1;
// resthcloseedgeenable = 1;
resthedgeanglefac = 1;
resthedgeangleenable = 0;
resthsurfmeshcurvfac = 1;
@ -1488,8 +1488,8 @@ void STLParameters :: Print (ostream & ost) const
<< ", fac = " << resthchartdistfac << endl
<< "line length: " << resthlinelengthenable
<< ", fac = " << resthlinelengthfac << endl
<< "close edges: " << resthcloseedgeenable
<< ", fac = " << resthcloseedgefac << endl
// << "close edges: " << resthcloseedgeenable
// << ", fac = " << resthcloseedgefac << endl
<< "edge angle: " << resthedgeangleenable
<< ", fac = " << resthedgeanglefac << endl;
}

View File

@ -282,8 +282,8 @@ public:
double resthchartdistfac = 1.2;
bool resthchartdistenable = true;
double resthcloseedgefac = 1.;
bool resthcloseedgeenable = true;
// double resthcloseedgefac = 1.;
// bool resthcloseedgeenable = true;
double resthedgeanglefac = 1.;
bool resthedgeangleenable = false;
@ -302,6 +302,13 @@ public:
void Print (ostream & ost) const;
};
inline ostream & operator<< (ostream & ost, const STLParameters & stlparam)
{
stlparam.Print (ost);
return ost;
}
void STLMeshing (STLGeometry & geom,
Mesh & mesh,

View File

@ -1029,6 +1029,9 @@ namespace netgen
else
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcol);
static Point<3> xa[129];
static Vec<3> na[129];
for (int hi = 0; hi < seia.Size(); hi++)
@ -1058,8 +1061,6 @@ namespace netgen
if (curv.IsHighOrder()) // && curv.IsSurfaceElementCurved(sei))
{
if (hoplotn > 128) hoplotn = 128;
Point<3> xa[129];
Vec<3> na[129];
for (int i = 0; i < hoplotn; i++)
{

View File

@ -29,6 +29,11 @@ namespace netgen
// vssolution.AddUserVisualizationObject (vis);
GetVSSolution().AddUserVisualizationObject (vis);
}
void DeleteUserVisualizationObject (UserVisualizationObject * vis)
{
// vssolution.AddUserVisualizationObject (vis);
GetVSSolution().DeleteUserVisualizationObject (vis);
}
VisualSceneSolution :: SolData :: SolData ()

View File

@ -233,7 +233,12 @@ public:
{
user_vis.Append (vis);
}
void DeleteUserVisualizationObject (UserVisualizationObject * vis)
{
int pos = user_vis.Pos(vis);
if (pos >= 0)
user_vis.Delete(pos);
}
private:
void GetClippingPlaneTrigs (NgArray<ClipPlaneTrig> & trigs, NgArray<ClipPlanePoint> & pts);

View File

@ -1131,7 +1131,7 @@ proc timer2 { } {
}
}
after 10 { timer2 }
after 40 { timer2 }
}
# after 1000 { timer2 }
timer2

View File

@ -537,6 +537,7 @@ namespace netgen
// delete ng_geometry;
// ng_geometry = hgeom;
ng_geometry = shared_ptr<NetgenGeometry> (hgeom);
geometryregister[i]->SetParameters(interp);
mesh.reset();
return TCL_OK;
@ -679,40 +680,96 @@ namespace netgen
int argc, tcl_const char *argv[])
{
char buf[20], lstring[200];
static int prev_np = -1;
static int prev_ne = -1;
static int prev_nse = -1;
if (mesh)
{
if (prev_np != mesh->GetNP())
{
sprintf (buf, "%u", unsigned(mesh->GetNP()));
Tcl_SetVar (interp, "::status_np", buf, 0);
prev_np = mesh->GetNP();
}
if (prev_ne != mesh->GetNE())
{
sprintf (buf, "%u", unsigned(mesh->GetNE()));
Tcl_SetVar (interp, "::status_ne", buf, 0);
prev_ne = mesh->GetNE();
}
if (prev_nse != mesh->GetNSE())
{
sprintf (buf, "%u", unsigned(mesh->GetNSE()));
Tcl_SetVar (interp, "::status_nse", buf, 0);
prev_nse = mesh->GetNSE();
}
auto tets_in_qualclass = mesh->GetQualityHistogram();
lstring[0] = 0;
for (int i = 0; i < tets_in_qualclass.Size(); i++)
{
sprintf (buf, " %d", tets_in_qualclass[i]);
strcat (lstring, buf);
}
for (int i = tets_in_qualclass.Size(); i < 20; i++)
strcat (lstring, " 0");
Tcl_SetVar (interp, "::status_tetqualclasses", lstring, 0);
}
else
{
if (prev_np != 0)
{
Tcl_SetVar (interp, "::status_np", "0", 0);
Tcl_SetVar (interp, "::status_ne", "0", 0);
Tcl_SetVar (interp, "::status_nse", "0", 0);
prev_np = 0;
}
if (prev_ne != 0)
{
Tcl_SetVar (interp, "::status_ne", "0", 0);
prev_ne = 0;
}
if (prev_nse != 0)
{
Tcl_SetVar (interp, "::status_nse", "0", 0);
prev_nse = 0;
}
Tcl_SetVar (interp, "::status_tetqualclasses", "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0", 0);
}
static string prev_working;
string working = multithread.running ? "working" : " ";
if (working != prev_working)
{
Tcl_SetVar (interp, "::status_working", working.c_str(), 0);
prev_working = working;
}
/*
if (multithread.running)
Tcl_SetVar (interp, "::status_working", "working", 0);
else
Tcl_SetVar (interp, "::status_working", " ", 0);
*/
Tcl_SetVar (interp, "::status_task", const_cast<char *>(multithread.task), 0);
sprintf (buf, "%lf", multithread.percent);
Tcl_SetVar (interp, "::status_percent", buf, 0);
lstring[0] = 0;
for (int i = 1; i <= tets_in_qualclass.Size(); i++)
static string prev_task;
if (prev_task != string(multithread.task))
{
sprintf (buf, " %d", tets_in_qualclass.Get(i));
strcat (lstring, buf);
prev_task = multithread.task;
Tcl_SetVar (interp, "::status_task", prev_task.c_str(), 0);
}
for (int i = tets_in_qualclass.Size()+1; i <= 20; i++)
strcat (lstring, " 0");
Tcl_SetVar (interp, "::status_tetqualclasses", lstring, 0);
static double prev_percent = -1;
if (prev_percent != multithread.percent)
{
prev_percent = multithread.percent;
sprintf (buf, "%lf", prev_percent);
Tcl_SetVar (interp, "::status_percent", buf, 0);
}
{
lock_guard<mutex> guard(tcl_todo_mutex);
@ -1226,6 +1283,10 @@ namespace netgen
mparam.parallel_meshing = atoi (Tcl_GetVar (interp, "::options.parallel_meshing", 0));
mparam.nthreads = atoi (Tcl_GetVar (interp, "::options.nthreads", 0));
if(atoi(Tcl_GetVar (interp, "::stloptions.resthcloseedgeenable", 0)))
mparam.closeedgefac = atof(Tcl_GetVar (interp, "::stloptions.resthcloseedgefac", 0));
else
mparam.closeedgefac = {};
//BaseMoveableMem::totalsize = 0;
// 1048576 * atoi (Tcl_GetVar (interp, "::options.memory", 0));

View File

@ -73,8 +73,6 @@ const char * ngscript[] = {""
,"set options.inverttets 0\n"
,"set options.inverttrigs 0\n"
,"set options.autozrefine 0\n"
,"set options.parallel_meshing 1\n"
,"set options.nthreads 4\n"
,"set options.meshsize 1000\n"
,"set options.minmeshsize 0\n"
,"set options.curvaturesafety 2\n"
@ -86,6 +84,8 @@ const char * ngscript[] = {""
,"set options.opterrpow 2\n"
,"set options.grading 0.5\n"
,"set options.printmsg 2\n"
,"set options.parallel_meshing 1\n"
,"set options.nthreads 4\n"
,"set debug.slowchecks 0\n"
,"set debug.debugoutput 0\n"
,"set debug.haltexistingline 0\n"
@ -1423,7 +1423,7 @@ const char * ngscript[] = {""
,"}\n"
,"}\n"
,"}\n"
,"after 10 { timer2 }\n"
,"after 40 { timer2 }\n"
,"}\n"
,"timer2\n"
,"proc bgerror { error } {\n"

View File

@ -7,7 +7,6 @@ if(WIN32)
$<TARGET_OBJECTS:interface>
$<TARGET_OBJECTS:geom2d>
$<TARGET_OBJECTS:csg>
$<TARGET_OBJECTS:stl>
$<TARGET_OBJECTS:visual>
$<TARGET_OBJECTS:occ>
@ -23,7 +22,7 @@ endif(WIN32)
add_library(nglib SHARED nglib.cpp ${nglib_objects})
if(NOT WIN32)
target_link_libraries( nglib PUBLIC mesh stl interface geom2d csg stl visual)
target_link_libraries( nglib PUBLIC mesh interface geom2d csg stl visual)
if(USE_GUI)
target_link_libraries( nglib PUBLIC stlvis geom2dvis csgvis )
endif(USE_GUI)

View File

@ -536,7 +536,7 @@ namespace nglib
Ng_Mesh * mesh,
int levels)
{
Refinement2d ref(*(SplineGeometry2d*)geom);
Refinement ref(*(SplineGeometry2d*)geom);
HPRefinement (*(Mesh*)mesh, &ref, levels);
}
@ -547,7 +547,7 @@ namespace nglib
Ng_Mesh * mesh,
int levels, double parameter)
{
Refinement2d ref(*(SplineGeometry2d*)geom);
Refinement ref(*(SplineGeometry2d*)geom);
HPRefinement (*(Mesh*)mesh, &ref, levels, parameter);
}
@ -856,8 +856,8 @@ namespace nglib
mp->Transfer_Parameters();
occparam.resthcloseedgeenable = mp->closeedgeenable;
occparam.resthcloseedgefac = mp->closeedgefact;
if(mp->closeedgeenable)
mparam.closeedgefac = mp->closeedgefact;
// Delete the mesh structures in order to start with a clean
// slate
@ -1090,7 +1090,7 @@ namespace nglib
// ------------------ Begin - Second Order Mesh generation functions ----------------
DLL_HEADER void Ng_Generate_SecondOrder(Ng_Mesh * mesh)
{
Refinement ref;
Refinement ref(*((Mesh*) mesh)->GetGeometry());
ref.MakeSecondOrder(*(Mesh*) mesh);
}
@ -1139,7 +1139,7 @@ namespace nglib
// ------------------ Begin - Uniform Mesh Refinement functions ---------------------
DLL_HEADER void Ng_Uniform_Refinement (Ng_Mesh * mesh)
{
Refinement ref;
Refinement ref(*((Mesh*)mesh)->GetGeometry());
ref.Refine ( * (Mesh*) mesh );
}

View File

@ -6,12 +6,10 @@ class _MeshsizeObject:
return MeshingParameters(curvaturesafety=1,
segmentsperedge=0.3,
grading=0.7,
surfcurvfac=0.25,
chartdistfac=0.8,
linelengthfac=0.2,
closeedgefac=0.5,
minedgelen=0.002,
edgeanglefac=0.25,
surfmeshcurvfac=1.,
optsteps3d=5)
@property
@ -19,12 +17,10 @@ class _MeshsizeObject:
return MeshingParameters(curvaturesafety=1.5,
segmentsperedge=0.5,
grading=0.5,
surfcurvfac=0.5,
chartdistfac=1,
linelengthfac=0.35,
closeedgefac=1,
minedgelen=0.02,
edgeanglefac=0.5,
surfmeshcurvfac=1.5,
optsteps3d=5)
@property
@ -32,12 +28,10 @@ class _MeshsizeObject:
return MeshingParameters(curvaturesafety=2,
segmentsperedge=1,
grading=0.3,
surfcurvfac=1.,
chartdistfac=1.5,
linelengthfac=0.5,
closeedgefac=2,
minedgelen=0.2,
edgeanglefac=1,
surfmeshcurvfac=2.,
optsteps3d=5)
@property
@ -45,12 +39,10 @@ class _MeshsizeObject:
return MeshingParameters(curvaturesafety=3,
segmentsperedge=2,
grading=0.2,
surfcurvfac=1.5,
chartdistfac=2,
linelengthfac=1.5,
closeedgefac=3.5,
minedgelen=1.,
edgeanglefac=1.5,
surfmeshcurvfac=3.,
optsteps3d=5)
@ -59,12 +51,10 @@ class _MeshsizeObject:
return MeshingParameters(curvaturesafety=5,
segmentsperedge=3,
grading=0.1,
surfcurvfac=3,
chartdistfac=5,
linelengthfac=3,
closeedgefac=5,
minedgelen=2.,
edgeanglefac=3.,
surfmeshcurvfac=5.,
optsteps3d=5)

View File

@ -1,5 +1,5 @@
FROM ubuntu:18.04
FROM ubuntu:19.10
ENV DEBIAN_FRONTEND=noninteractive
MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk clang-tidy python3-distutils clang liboce-ocaf-dev
RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk clang-tidy python3-distutils clang libocct-data-exchange-dev
ADD . /root/src/netgen

View File

@ -0,0 +1,91 @@
import json
import sys
import subprocess
import statistics
def readData(a, files):
amin=[]
amax=[]
amin1=[]
amax1=[]
bad=[]
ne1d=[]
ne2d=[]
ne3d=[]
for f in files:
for t in a[f]:
if t['ne1d']>0:
ne1d.append(t['ne1d'])
if t['ne2d']>0:
ne2d.append(t['ne2d'])
if t['ne3d']>0:
ne3d.append(t['ne3d'])
if t['total_badness']>0.0:
bad.append(t['total_badness'])
if 'angles_tet' in t:
amin.append(t['angles_tet'][0])
amax.append(t['angles_tet'][1])
if 'angles_trig' in t:
amin1.append(t['angles_trig'][0])
amax1.append(t['angles_trig'][1])
return {
"min tet angle":amin,
"max tet angle" : amax,
"min trig angle":amin1,
"max trig angle" : amax1,
"badness" : bad,
"#edges" : ne1d,
"#trigs" : ne2d,
"#tets" : ne3d,
}
import matplotlib.pyplot as plt
ref = 'master'
if len(sys.argv)>1:
ref = sys.argv[1]
res = subprocess.run(['git','show','{}:./results.json'.format(ref)], capture_output=True)
s = json.loads(res.stdout.decode())
if len(sys.argv) > 2:
ref2 = sys.argv[2]
res = subprocess.run(['git','show','{}:./results.json'.format(ref2)], capture_output=True)
s2 = res.stdout.decode()
else:
ref2 = 'current'
s2 = open('results.json','r').read()
s2 = json.loads(s2)
filenames = [f for f in s if f in s2]
data = readData(s, filenames)
data2 = readData(s2, filenames)
n = len(data)+1
fig,ax = plt.subplots(figsize=(10,7))
for i,d in enumerate(['min trig angle','min tet angle','max trig angle','max tet angle']):
ax = plt.subplot(2,5,i+1)
plt.title(d)
ax.set_xticks([1,2])
if len(data[d])==0 or len(data2[d])==0:
continue
plt.violinplot([data[d],data2[d]], showmedians=True)
med = statistics.median(data[d])
plt.hlines(med, 1,2, linestyle='dotted')
if d=='badness':
ax.set_yscale('log')
ax.set_xticklabels([ref, ref2])
for i,d in enumerate(['badness','#edges','#trigs','#tets']):
ax = plt.subplot(2,5,6+i)
plt.title('difference '+d+' (in %)')
# plt.violinplot([(y-x)/x for x,y in zip(data[d],data2[d])], showmedians=True)
plt.boxplot([100.0*(y-x)/x for x,y in zip(data[d],data2[d])])
plt.hlines(0.0, 0.5,1.5, linestyle='dotted')
# plt.savefig('comparison.png', dpi=100)
plt.show()

File diff suppressed because it is too large Load Diff

3442
tests/pytest/results.json Normal file

File diff suppressed because it is too large Load Diff

View File

@ -1,36 +0,0 @@
number_elements = {}
number_elements['cylsphere.geo'] = (706,231,490,706,2797,17554)
number_elements['cubeandspheres.geo'] = (98,100,98,98,366,1078)
number_elements['ellipsoid.geo'] = (1271,551,595,1268,5607,38173)
number_elements['manyholes2.geo'] = (128244)
number_elements['sculpture.geo'] = (477,140,260,477,1330,6769)
number_elements['ortho.geo'] = (6,6,6,6,34,180)
number_elements['ellipticcone.geo'] = (4973,573,1765,4918,13410,70483)
number_elements['cube.geo'] = (6,6,6,6,28,178)
number_elements['twobricks.geo'] = (42,22,22,42,177,595)
number_elements['revolution.geo'] = (8310,1249,3856,8269,33078,202941)
number_elements['circle_on_cube.geo'] = (636,39,189,631,2035,12237)
number_elements['sphereincube.geo'] = (508,173,339,515,1652,13829)
number_elements['twocubes.geo'] = (42,22,22,42,177,595)
number_elements['boundarycondition.geo'] = (39,22,22,39,165,508)
number_elements['ellipticcyl.geo'] = (2202,324,1106,2190,8245,55199)
number_elements['trafo.geo'] = (5154,1358,2389,5141,17948,92850)
number_elements['boxcyl.geo'] = (843,146,364,843,3700,18677)
number_elements['sphere.geo'] = (126,56,80,126,347,2342)
number_elements['torus.geo'] = (5520,2171,2739,5510,25402,177967)
number_elements['shaft.geo'] = (2609,808,1666,2594,11226,64172)
number_elements['cone.geo'] = (1215,447,678,1211,4404,27336)
number_elements['cubeandring.geo'] = (2014,231,612,1988,7671,38095)
number_elements['manyholes.geo'] = (176503,28896,70408)
number_elements['period.geo'] = (3263,574,1349,3236,11645,68523)
number_elements['lshape3d.geo'] = (18,12,12,18,93,335)
number_elements['cubemsphere.geo'] = (4708,773,1460,4667,17655,114554)
number_elements['twocyl.geo'] = (578,147,403,578,1887,13537)
number_elements['cubemcyl.geo'] = (19712,3225,8153,19438,89202,524684)
number_elements['matrix.geo'] = (5207,1888,2790,5149,16205,101146)
number_elements['fichera.geo'] = (35,18,18,35,209,496)
number_elements['cylinder.geo'] = (404,101,256,404,1161,8076)
number_elements['part1.stl'] = (1228,620,727,1216,1548,3498)
number_elements['hinge.stl'] = (1995,1399,1532,1987,2881,4621)
number_elements['frame.step'] = (195213,30333,58955)
number_elements['screw.step'] = (2021,7011,23730)

View File

@ -3,36 +3,64 @@ import os, pytest
from netgen.meshing import meshsize, MeshingParameters, SetMessageImportance
import netgen.csg as csg
import netgen.stl as stl
import netgen.geom2d as geom2d
from pyngcore import TaskManager
import json
try:
import netgen.occ as occ
has_occ = True
has_occ = occ.occ_version >= "7.4.0"
except ImportError:
has_occ = False
from results import *
SetMessageImportance(0)
def round(x, digits=11):
try:
return float(("{:."+str(digits)+"g}").format(x))
except: #list
return [float(("{:."+str(digits)+"g}").format(y)) for y in x]
def getData(mesh, mp):
out = {}
out['ne1d'] = len(mesh.Elements1D())
out['ne2d'] = len(mesh.Elements2D())
out['ne3d'] = len(mesh.Elements3D())
# round badness to avoid fluctuations in last digits
out["total_badness"] = round(mesh.CalcTotalBadness(mp))
angles = mesh.CalcMinMaxAngle()
out["angles_trig"] = round(angles["trig"], 5)
out["angles_tet"] = round(angles["tet"], 5)
out["quality_histogram"] = str(list(mesh.GetQualityHistogram()))
return out
def checkData(mesh, mp, ref):
data = getData(mesh, mp)
assert ref['ne1d'] == data['ne1d']
assert ref['ne2d'] == data['ne2d']
assert ref['ne3d'] == data['ne3d']
assert json.loads(ref['quality_histogram']) == pytest.approx(json.loads(data['quality_histogram']), abs=1, rel=0.4)
assert ref['total_badness'] == pytest.approx(data['total_badness'], rel=1e-5)
assert ref['angles_trig'] == pytest.approx(data['angles_trig'], rel=1e-4)
assert ref['angles_tet'] == pytest.approx(data['angles_tet'], rel=1e-4)
# get tutorials
def getFiles(fileEnding):
r, d, files = next(os.walk(os.path.join("..","..","tutorials")))
return (f for f in files if f.endswith(fileEnding))
return [f for f in files if f.endswith(fileEnding)]
# get additional tests
def getAdditionalFiles(fileEnding):
r, d, files = next(os.walk("geofiles"))
return [f for f in files if f.endswith(fileEnding)]
def getCheckFunc(filename):
def func(mesh,i):
if filename in number_elements:
# number of elements should be in 2% range of expected value
assert mesh.ne == pytest.approx(number_elements[filename][i], rel=0.02)
return func
@pytest.fixture
def refdata():
return json.load(open('results.json','r'))
def getResultFunc(filename):
def resultFunc(mesh):
results = {}
results["number_elements"] = mesh.ne
return results
return resultFunc
def getMeshingparameters(filename):
standard = [{}] + [{ "mp" : ms } for ms in (meshsize.very_coarse, meshsize.coarse, meshsize.moderate, meshsize.fine, meshsize.very_fine)]
standard = [MeshingParameters()] + [MeshingParameters(ms) for ms in (meshsize.very_coarse, meshsize.coarse, meshsize.moderate, meshsize.fine, meshsize.very_fine)]
if filename == "shell.geo":
return [] # do not test this example cause it needs so long...
if filename == "extrusion.geo":
@ -43,61 +71,85 @@ def getMeshingparameters(filename):
return standard[:3] # this gets too big for finer meshsizes
if filename == "screw.step":
return standard[3:] # coarser meshes don't work here
if filename == "cylsphere.geo":
return standard[0:2] + standard[3:] # coarse gives inconsistent reults (other mesh on MacOS)
if filename == "part1.stl":
return standard[0:1] + standard[2:] # very coarse does not work
return standard
_geofiles = [f for f in getFiles(".geo")] + [f for f in getFiles(".stl")]
_geofiles = getFiles(".in2d") + getFiles(".geo") + getFiles(".stl")
if has_occ:
_geofiles += [f for f in getFiles(".step")]
_geofiles += getFiles(".step")
_geofiles.sort()
_additional_testfiles = getAdditionalFiles(".stl")
if has_occ:
_additional_testfiles += getAdditionalFiles(".step")
_additional_testfiles.sort()
def generateMesh(filename, mp):
folder = os.path.join("..","..","tutorials") if filename in _geofiles else "geofiles"
if filename.endswith(".geo"):
geo = csg.CSGeometry(os.path.join("..","..","tutorials", filename))
geo = csg.CSGeometry(os.path.join(folder, filename))
elif filename.endswith(".stl"):
geo = stl.STLGeometry(os.path.join("..","..","tutorials", filename))
geo = stl.STLGeometry(os.path.join(folder, filename))
elif filename.endswith(".step"):
geo = occ.OCCGeometry(os.path.join("..","..","tutorials", filename))
return geo.GenerateMesh(**mp)
geo = occ.OCCGeometry(os.path.join(folder, filename))
elif filename.endswith(".in2d"):
geo = geom2d.SplineGeometry(os.path.join(folder, filename))
return geo.GenerateMesh(mp)
def isSlowTest(filename):
return filename in ["cubemcyl.geo", "frame.step", "revolution.geo", "manyholes.geo", "torus.geo",
"cubemsphere.geo", "manyholes2.geo", "matrix.geo", "trafo.geo", "ellipticcone.geo",
"period.geo", "shaft.geo", "cubeandring.geo", "ellipticcyl.geo",
"ellipsoid.geo", "cone.geo"]
"ellipsoid.geo", "cone.geo", "plane.stl"]
def getParamForTest(filename):
return pytest.param(filename, getCheckFunc(filename), marks=pytest.mark.slow) if isSlowTest(filename) \
else (filename, getCheckFunc(filename))
def getParameters():
res = []
for f in _geofiles + _additional_testfiles:
for i,mp in enumerate(getMeshingparameters(f)):
if isSlowTest(f):
res.append( pytest.param(f, mp, i, marks=pytest.mark.slow ) )
else:
res.append( (f, mp, i) )
return res
@pytest.mark.parametrize(("filename, checkFunc"), [getParamForTest(f) for f in _geofiles])
def test_geoFiles(filename, checkFunc):
import filecmp, pyngcore
for i, mp in enumerate(getMeshingparameters(filename)):
@pytest.mark.parametrize(("filename", "mp", "i"), getParameters())
def test_geoFiles(filename, mp, i, refdata):
ref = refdata[filename]
import filecmp
print("load geo", filename)
mp = MeshingParameters(mp, parallel_meshing=False)
mesh = generateMesh(filename, mp)
if checkFunc is not None:
checkFunc(mesh,i)
mesh.Save(filename+'_seq.vol.gz')
with pyngcore.TaskManager():
with TaskManager():
mesh_par = generateMesh(filename, mp)
mesh_par.Save(filename+'_par.vol.gz')
assert filecmp.cmp(filename+'_seq.vol.gz', filename+'_par.vol.gz')
checkData(mesh, mp, ref[i])
import time
def generateResultFile():
with open("results.py", "w") as f:
print("number_elements = {}", file=f)
for _file, _func in ((gf, getResultFunc(gf)) for gf in _geofiles):
import re, time
data = {}
with TaskManager():
for _file in _geofiles + _additional_testfiles:
print("generate "+_file)
start = time.time()
print("write", _file)
mps = getMeshingparameters(_file)
if not mps:
continue
results = [_func(generateMesh(_file, mp)) for mp in mps]
print("number_elements['{}'] = {}".format(_file, "(" + ",".join((str(r["number_elements"]) for r in results)) + ")"), file=f)
meshdata = []
for mp in mps:
mesh = generateMesh(_file, mp)
meshdata.append( getData(mesh, mp) )
data[_file] = meshdata
print("needed", time.time() - start, "seconds")
s = json.dumps(data, sort_keys=True, indent=4)
open("results.json", "w").write(s)
print("done")
if __name__ == "__main__":
generateResultFile()