ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-05-18 16:20:50 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' of https://github.com/NGSolve/netgen

Browse Source
This commit is contained in:
Duncan McDougall 2020-06-26 11:10:10 +01:00
commit b4337c5df9
29 changed files with 860 additions and 160 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
CMakeLists.txt
View File

@ -35,8 +35,6 @@ option( USE_SUPERBUILD "use ccache" ON)
set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_modules")
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
if(APPLE)
set(INSTALL_DIR_DEFAULT /Applications/Netgen.app)
else(APPLE)
@ -78,6 +76,8 @@ else()
endif()
endif()
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
include (${CMAKE_CURRENT_LIST_DIR}/cmake/generate_version_file.cmake)
set(CPACK_PACKAGE_VERSION "${NETGEN_VERSION}")
@ -191,7 +191,6 @@ check_include_files (dlfcn.h HAVE_DLFCN_H)
if(HAVE_DLFCN_H)
add_definitions(-DHAVE_DLFCN_H)
endif()
add_definitions(-DNETGEN_VERSION="${NETGEN_VERSION}")
include_directories(BEFORE ${CMAKE_CURRENT_BINARY_DIR})

2
cmake/cmake_modules/FindOpenCasCade.cmake
View File

@ -27,7 +27,7 @@ endif(WIN32)
if(OCC_LIBRARY AND NOT OCC_LIBRARY_DIR)
get_filename_component(OCC_LIBRARY_DIR ${OCC_LIBRARY} PATH)
endif(OCC_LIBRARY)
endif(OCC_LIBRARY AND NOT OCC_LIBRARY_DIR)
if(OCC_INCLUDE_DIR)
file(STRINGS ${OCC_INCLUDE_DIR}/Standard_Version.hxx OCC_MAJOR

8
libsrc/core/array.hpp
View File

@ -1284,7 +1284,7 @@ namespace ngcore
/// bubble sort array
template <class T, class S>
inline void BubbleSort (FlatArray<T> data, FlatArray<S> slave)
inline void BubbleSort (FlatArray<T> data, FlatArray<S> index)
{
for (size_t i = 0; i < data.Size(); i++)
for (size_t j = i+1; j < data.Size(); j++)
@ -1294,9 +1294,9 @@ namespace ngcore
data[i] = data[j];
data[j] = hv;
S hvs = slave[i];
slave[i] = slave[j];
slave[j] = hvs;
S hvs = index[i];
index[i] = index[j];
index[j] = hvs;
}
}

1
libsrc/csg/csgeom.cpp
View File

@ -130,6 +130,7 @@ namespace netgen
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)
{

20
libsrc/general/ngarray.hpp
View File

@ -730,7 +730,7 @@ namespace netgen
/// bubble sort array
template <class T, class S>
inline void BubbleSort (NgFlatArray<T> & data, NgFlatArray<S> & slave)
inline void BubbleSort (NgFlatArray<T> & data, NgFlatArray<S> & index)
{
for (int i = 0; i < data.Size(); i++)
for (int j = i+1; j < data.Size(); j++)
@ -740,16 +740,16 @@ namespace netgen
data[i] = data[j];
data[j] = hv;
S hvs = slave[i];
slave[i] = slave[j];
slave[j] = hvs;
S hvs = index[i];
index[i] = index[j];
index[j] = hvs;
}
}
template <class T, class S>
void QuickSortRec (NgFlatArray<T> & data,
NgFlatArray<S> & slave,
NgFlatArray<S> & index,
int left, int right)
{
int i = left;
@ -764,20 +764,20 @@ namespace netgen
if (i <= j)
{
ngcore::Swap (data[i], data[j]);
ngcore::Swap (slave[i], slave[j]);
ngcore::Swap (index[i], index[j]);
i++; j--;
}
}
while (i <= j);
if (left < j) QuickSortRec (data, slave, left, j);
if (i < right) QuickSortRec (data, slave, i, right);
if (left < j) QuickSortRec (data, index, left, j);
if (i < right) QuickSortRec (data, index, i, right);
}
template <class T, class S>
void QuickSort (NgFlatArray<T> & data, NgFlatArray<S> & slave)
void QuickSort (NgFlatArray<T> & data, NgFlatArray<S> & index)
{
if (data.Size() > 1)
QuickSortRec (data, slave, 0, data.Size()-1);
QuickSortRec (data, index, 0, data.Size()-1);
}

13
libsrc/geom2d/genmesh2d.cpp
View File

@ -527,11 +527,16 @@ namespace netgen
for (PointIndex pix = nextpi[c1], ix = 0; pix != c2; pix = nextpi[pix], ix++)
{
Point<3> px = (*mesh)[pix];
for (PointIndex piy = nextpi[c2], iy = 0; piy != c3; piy = nextpi[piy], iy++)
{
Point<3> p = (*mesh)[pix] + ( (*mesh)[piy] - (*mesh)[c2] );
pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (p , 1, FIXEDPOINT);
double lam = Dist((*mesh)[piy],(*mesh)[c2]) / Dist((*mesh)[c3],(*mesh)[c2]);
auto pix1 = pts[(nex+1)*ney+ix+1];
auto pnew = px + lam*((*mesh)[pix1]-px);
pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (pnew, 1, FIXEDPOINT);
}
}
for (int i = 0; i < ney; i++)
for (int j = 0; j < nex; j++)
@ -545,6 +550,10 @@ namespace netgen
mesh -> AddSurfaceElement (el);
}
char* material;
geometry.GetMaterial(domnr, material);
if(material)
mesh->SetMaterial(domnr, material);
}

15
libsrc/geom2d/geometry2d.hpp
View File

@ -133,7 +133,7 @@ namespace netgen
NgArray<char*> materials;
NgArray<double> maxh;
NgArray<bool> quadmeshing;
NgArray<bool> tensormeshing;
Array<bool> tensormeshing;
NgArray<int> layer;
NgArray<string*> bcnames;
double elto0 = 1.0;
@ -216,9 +216,20 @@ namespace netgen
}
bool GetDomainTensorMeshing ( int domnr )
{
if ( tensormeshing.Size() ) return tensormeshing[domnr-1];
if ( tensormeshing.Size()>=domnr ) return tensormeshing[domnr-1];
else return false;
}
void SetDomainTensorMeshing ( int domnr, bool tm )
{
if ( tensormeshing.Size()<domnr )
{
auto oldsize = tensormeshing.Size();
tensormeshing.SetSize(domnr);
for(auto i : IntRange(oldsize, domnr-1))
tensormeshing[i] = false;
}
tensormeshing[domnr-1] = tm;
}
int GetDomainLayer ( int domnr )
{
if ( layer.Size() ) return layer[domnr-1];

1
libsrc/geom2d/python_geom2d.cpp
View File

@ -392,6 +392,7 @@ DLL_HEADER void ExportGeom2d(py::module &m)
}, py::arg("mp") = nullptr,
py::call_guard<py::gil_scoped_release>(),
meshingparameter_description.c_str())
.def("_SetDomainTensorMeshing", &SplineGeometry2d::SetDomainTensorMeshing)
;
}

8
libsrc/include/nginterface_v2_impl.hpp
View File

@ -111,7 +111,13 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const
ret.faces.num = 0;
ret.faces.ptr = NULL;
if (mesh->GetDimension() == 2)
if (mesh->GetDimension() == 3)
{
ret.facets.num = 0;
ret.facets.base = 0;
ret.facets.ptr = nullptr;
}
else if (mesh->GetDimension() == 2)
{
ret.facets.num = 1;
ret.facets.base = 0;

18
libsrc/interface/readtetmesh.cpp
View File

@ -154,7 +154,7 @@ namespace netgen
break;
case 7:
// NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), PID:
// NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), PID:
{
cout << "read nodes" << endl;
for(int i=0; i<nnodes; i++)
@ -175,7 +175,7 @@ namespace netgen
break;
case 9:
// MasterNodeID, SlaveNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2)
// MasterNodeID, MinionNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2)
for(int i=0; i<nperiodicmasternodes; i++)
{
for(int j=0; j<2; j++)
@ -191,7 +191,7 @@ namespace netgen
break;
case 11:
// MasterNodeID, 3-SlaveNodeID's, 3-TranslCodes (1=dS1 2=dS2 3=dS1+dS2)
// MasterNodeID, 3-MinionNodeID's, 3-TranslCodes (1=dS1 2=dS2 3=dS1+dS2)
for(int i=0; i<ncornerperiodicmasternodes; i++)
{
for(int j=0; j<4; j++)
@ -208,7 +208,7 @@ namespace netgen
break;
case 13:
//MasterNodeID, 7-SlaveNodeID's, TranslCodes
//MasterNodeID, 7-MinionNodeID's, TranslCodes
for(int i=0; i<ncubicperiodicmasternodes; i++)
{
for(int j=0; j<8; j++)
@ -220,7 +220,7 @@ namespace netgen
break;
case 14:
// EdgeID, NodeID0, NodeID1, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), PID
// EdgeID, NodeID0, NodeID1, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), PID
cout << "read edges" << endl;
// nullstarted = false;
segmentdata.SetSize(nedges);
@ -243,7 +243,7 @@ namespace netgen
break;
case 16:
// MasterEdgeID, SlaveEdgeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2)
// MasterEdgeID, MinionEdgeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2)
for(int i=0; i<nperiodicmasteredges; i++)
in >> dummyint >> dummyint >> dummyint;
break;
@ -254,7 +254,7 @@ namespace netgen
break;
case 18:
// MasterEdgeID, 3 SlaveEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2)
// MasterEdgeID, 3 MinionEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2)
for(int i=0; i<ncornerperiodicmasteredges; i++)
{
in >> dummyint;
@ -266,7 +266,7 @@ namespace netgen
break;
case 19:
// FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PSlave), PID
// FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PMinion), PID
{
//Segment seg;
int segnum_ng[3];
@ -343,7 +343,7 @@ namespace netgen
break;
case 21:
// MasterFaceID, SlaveFaceID, TranslCode (1=dS1 2=dS2)
// MasterFaceID, MinionFaceID, TranslCode (1=dS1 2=dS2)
{
Vec<3> randomvec(-1.32834,3.82399,0.5429151);
int maxtransl = -1;

21
libsrc/interface/writeabaqus.cpp
View File

@ -268,25 +268,25 @@ void WriteAbaqusFormat (const Mesh & mesh,
cout << "masternode = " << masternode << " = "
<< mesh.Point(masternode) << endl;
NgArray<int> slaves(3);
NgArray<int> minions(3);
for (i = 1; i <= 3; i++)
{
mesh.GetIdentifications().GetPairs (i, pairs);
for (j = 1; j <= pairs.Size(); j++)
{
if (pairs.Get(j).I1() == masternode)
slaves.Elem(i) = pairs.Get(j).I2();
minions.Elem(i) = pairs.Get(j).I2();
}
cout << "slave(" << i << ") = " << slaves.Get(i)
<< " = " << mesh.Point(slaves.Get(i)) << endl;
cout << "minion(" << i << ") = " << minions.Get(i)
<< " = " << mesh.Point(minions.Get(i)) << endl;
}
outfile << "**\n"
<< "*NSET,NSET=CTENODS\n"
<< slaves.Get(1) << ", "
<< slaves.Get(2) << ", "
<< slaves.Get(3) << endl;
<< minions.Get(1) << ", "
<< minions.Get(2) << ", "
<< minions.Get(3) << endl;
outfile << "**\n"
@ -300,7 +300,7 @@ void WriteAbaqusFormat (const Mesh & mesh,
<< "*BOUNDARY, OP=NEW\n";
for (j = 1; j <= 3; j++)
{
Vec3d v(mesh.Point(masternode), mesh.Point(slaves.Get(j)));
Vec3d v(mesh.Point(masternode), mesh.Point(minions.Get(j)));
double vlen = v.Length();
int dir = 0;
if (fabs (v.X()) > 0.9 * vlen) dir = 2;
@ -308,7 +308,7 @@ void WriteAbaqusFormat (const Mesh & mesh,
if (fabs (v.Z()) > 0.9 * vlen) dir = 1;
if (!dir)
cout << "ERROR: Problem with rigid body constraints" << endl;
outfile << slaves.Get(j) << ", " << dir << ",, 0.\n";
outfile << minions.Get(j) << ", " << dir << ",, 0.\n";
}
outfile << "**\n"
@ -333,14 +333,13 @@ void WriteAbaqusFormat (const Mesh & mesh,
mpc << "4" << "\n";
mpc << pairs.Get(j).I2() << "," << k << ", -1.0, ";
mpc << pairs.Get(j).I1() << "," << k << ", 1.0, ";
mpc << slaves.Get(i) << "," << k << ", 1.0, ";
mpc << minions.Get(i) << "," << k << ", 1.0, ";
mpc << masternode << "," << k << ", -1.0 \n";
}
}
}
}
cout << "done" << endl;
}

8
libsrc/interface/writefeap.cpp
View File

@ -144,11 +144,11 @@ void WriteFEAPFormat (const Mesh & mesh,
// BEGIN CONTACT OUTPUT
/*
int masterindex, slaveindex;
int masterindex, minionindex;
cout << "Master Surface index = ";
cin >> masterindex;
cout << "Slave Surface index = ";
cin >> slaveindex;
cout << "Minion Surface index = ";
cin >> minionindex;
// CONTACT SURFACE 1
@ -196,7 +196,7 @@ void WriteFEAPFormat (const Mesh & mesh,
Element2d sel = mesh.SurfaceElement(i);
if (invertsurf)
sel.Invert();
if (mesh.GetFaceDescriptor(sel.GetIndex ()).BCProperty() == slaveindex)
if (mesh.GetFaceDescriptor(sel.GetIndex ()).BCProperty() == minionindex)
{
zz++;
outfile.width(14);

18
libsrc/interface/writetet.cpp
View File

@ -427,7 +427,7 @@ namespace netgen
<< numedges << " "
<< numnodes << endl << endl;
outfile << "// NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), "<< uidpid <<":\n" \
outfile << "// NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), "<< uidpid <<":\n" \
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
@ -515,7 +515,7 @@ namespace netgen
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n" \
<< n2 << "\n" \
<< "\n" \
<< "// MasterNodeID, SlaveNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n" \
<< "// MasterNodeID, MinionNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n" \
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
for(int i=0; i<id_groups.Size(); i++)
{
@ -538,7 +538,7 @@ namespace netgen
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n" \
<< n4 << "\n" \
<< "\n" \
<< "// MasterNodeID, 3-SlaveNodeID's, 3-TranslCodes (1=dS1 2=dS2 3=dS1+dS2):\n" \
<< "// MasterNodeID, 3-MinionNodeID's, 3-TranslCodes (1=dS1 2=dS2 3=dS1+dS2):\n" \
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
@ -565,7 +565,7 @@ namespace netgen
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n" \
<< n8 << "\n" \
<< "\n" \
<< "// MasterNodeID, 7-SlaveNodeID's, TranslCodes:\n" \
<< "// MasterNodeID, 7-MinionNodeID's, TranslCodes:\n" \
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
for(int i=0; i<id_groups.Size(); i++)
{
@ -586,7 +586,7 @@ namespace netgen
outfile << "// EdgeID, NodeID0, NodeID1, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), "<<uidpid<<":\n" \
outfile << "// EdgeID, NodeID0, NodeID1, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), "<<uidpid<<":\n" \
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
@ -760,7 +760,7 @@ namespace netgen
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"\
<< n2 << "\n" \
<< "\n"\
<< "// MasterEdgeID, SlaveEdgeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n"\
<< "// MasterEdgeID, MinionEdgeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n"\
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
for(int i=0; i<id_groups.Size(); i++)
{
@ -782,7 +782,7 @@ namespace netgen
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"\
<< n4 << "\n" \
<< "\n"\
<< "// MasterEdgeID, 3 SlaveEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n"\
<< "// MasterEdgeID, 3 MinionEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n"\
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
for(int i=0; i<id_groups.Size(); i++)
{
@ -801,7 +801,7 @@ namespace netgen
outfile << endl;
outfile << "// FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PSlave), "<<uidpid<<":\n" \
outfile << "// FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PMinion), "<<uidpid<<":\n" \
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
@ -921,7 +921,7 @@ namespace netgen
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"\
<< n2 << "\n" \
<< "\n"\
<< "// MasterFaceID, SlaveFaceID, TranslCode (1=dS1 2=dS2):\n"\
<< "// MasterFaceID, MinionFaceID, TranslCode (1=dS1 2=dS2):\n"\
<< "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
for(int i=0; i<id_groups.Size(); i++)
{

4
libsrc/meshing/CMakeLists.txt
View File

@ -12,7 +12,7 @@ add_library(mesh ${NG_LIB_TYPE}
smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp
topology.cpp triarls.cpp validate.cpp bcfunctions.cpp
parallelmesh.cpp paralleltop.cpp paralleltop.hpp basegeom.cpp
python_mesh.cpp hexarls.cpp
python_mesh.cpp hexarls.cpp surfacegeom.cpp
../../ng/onetcl.cpp
${mesh_object_libs}
)
@ -37,6 +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
python_mesh.hpp surfacegeom.hpp
DESTINATION ${NG_INSTALL_DIR_INCLUDE}/meshing COMPONENT netgen_devel
)

2
libsrc/meshing/bcfunctions.cpp
View File

@ -322,7 +322,7 @@ namespace netgen
colours_sorted.SetSize(all_colours.Size()+1);
faces_sorted = 0;
// Slave NgArray to identify the colours the faces were assigned to,
// Index NgArray to identify the colours the faces were assigned to,
// after the bubble sort routine to sort the automatic boundary
// identifiers according to the number of surface mesh elements
// of a given colour

2
libsrc/meshing/bcfunctions.hpp
View File

@ -26,7 +26,7 @@ namespace netgen
- Use colour index 0 (zero) for all faces with no colour defined
- Calculate the number of faces of the surface mesh for each colour
- Sort the number of surface elements in ascending order, with the
colour indices as a slave
colour indices as a index
- Use the indices of the sorted array as the BC property number
Example: If there are 3 colours, present in the mesh and the number

134
libsrc/meshing/boundarylayer.cpp
View File

@ -159,10 +159,13 @@ namespace netgen
auto& old_fd = mesh.GetFaceDescriptor(si);
int domout = blp.outside ? old_fd.DomainOut() : blp.new_matnrs[layer-1];
int domin = blp.outside ? blp.new_matnrs[layer-1] : old_fd.DomainIn();
FaceDescriptor fd(max_surface_index-1,
// -1 surf nr is so that curving does not do anything
FaceDescriptor fd(-1,
domin, domout, -1);
fd.SetBCProperty(max_surface_index);
mesh.AddFaceDescriptor(fd);
mesh.SetBCName(max_surface_index-1,
"mapped_" + old_fd.GetBCName());
return max_surface_index;
}
return last_layer_surface_index_map[si];
@ -269,55 +272,53 @@ namespace netgen
if(blp.grow_edges)
for(SegmentIndex sei = 0; sei < nseg; sei++)
{
PointIndex seg_p1 = mesh[sei][0];
PointIndex seg_p2 = mesh[sei][1];
auto& segi = mesh[sei];
// Only go in if the segment is still active, and if both its
// surface index is part of the "hit-list"
if(segsel.Test(sei))
{
if(blp.surfid.Contains(mesh[sei].si))
{
// clear the bit to indicate that this segment has been processed
segsel.Clear(sei);
// Find matching segment pair on other surface
for(SegmentIndex sej = 0; sej < nseg; sej++)
if(blp.surfid.Contains(segi.si))
{
PointIndex segpair_p1 = mesh[sej][1];
PointIndex segpair_p2 = mesh[sej][0];
// clear the bit to indicate that this segment has been processed
segsel.Clear(sei);
// Find the segment pair on the neighbouring surface element
// Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0]
if(segsel.Test(sej) && ((segpair_p1 == seg_p1) && (segpair_p2 == seg_p2)))
// Find matching segment pair on other surface
for(SegmentIndex sej = 0; sej < nseg; sej++)
{
// clear bit to indicate that processing of this segment is done
segsel.Clear(sej);
auto& segj = mesh[sej];
// Find the segment pair on the neighbouring surface element
// Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0]
if(segsel.Test(sej) && ((segi[0] == segj[1]) && (segi[1] == segj[0])))
{
// clear bit to indicate that processing of this segment is done
segsel.Clear(sej);
// Only worry about those surfaces which are not in the
// boundary layer list
if(!blp.surfid.Contains(mesh[sej].si))
// if segj is not in surfel list we nned to add quads
if(!blp.surfid.Contains(segj.si))
{
SurfaceElementIndex pnt_commelem;
SetInvalid(pnt_commelem);
auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segpair_p1);
auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2);
auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segj[0]);
auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segj[1]);
for(auto pnt1_sei : pnt1_elems)
if(mesh[pnt1_sei].GetIndex() == mesh[sej].si)
if(mesh[pnt1_sei].GetIndex() == segj.si)
for(auto pnt2_sei : pnt2_elems)
if(pnt1_sei == pnt2_sei)
pnt_commelem = pnt1_sei;
if(IsInvalid(pnt_commelem))
throw Exception("Couldn't find element on other side for " + ToString(segpair_p1) + " to " + ToString(segpair_p2));
throw Exception("Couldn't find element on other side for " + ToString(segj[0]) + " to " + ToString(segj[1]));
const auto& commsel = mesh[pnt_commelem];
auto surfelem_vect = GetSurfaceNormal(mesh, commsel);
if(blp.outside)
surfelem_vect *= -1;
Element2d sel(QUAD);
auto seg_p1 = segi[0];
auto seg_p2 = segi[1];
if(blp.outside)
Swap(seg_p1, seg_p2);
sel[0] = seg_p1;
@ -330,7 +331,7 @@ namespace netgen
{
domains_to_surf_index[domains] = ++max_surface_index;
domains_to_surf_index[make_tuple(max_surface_index, get<1>(domains), get<2>(domains))] = max_surface_index;
FaceDescriptor fd(max_surface_index-1,
FaceDescriptor fd(-1,
get<1>(domains),
get<2>(domains),
-1);
@ -365,46 +366,37 @@ namespace netgen
seg_2.edgenr = pi_to_edgenr[points];
seg_2.si = new_index;
mesh.AddSegment(seg_2);
mesh[sej].si = new_index;
}
// in last layer insert new segments
if(layer == blp.heights.Size())
{
Segment s1 = mesh[sei];
Segment s2 = mesh[sej];
Segment s1 = segi;
Segment s2 = segj;
s1.edgenr = ++max_edge_nr;
s2.edgenr = max_edge_nr;
bool create_it = true;
if(blp.surfid.Contains(mesh[sej].si))
{
if(last_layer_surface_index_map.find(s1.si) != last_layer_surface_index_map.end() &&
last_layer_surface_index_map.find(s2.si) != last_layer_surface_index_map.end())
// edge already mapped
create_it = false;
s2.si = map_surface_index(s2.si);
}
if(blp.surfid.Contains(segj.si))
s2.si = map_surface_index(segj.si);
else
{
s2.si = domains_to_surf_index[make_tuple(s2.si,
blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())];
auto side_surf = domains_to_surf_index[make_tuple(s2.si, blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())];
if(blp.outside)
s2.si = side_surf;
else
segj.si = side_surf;
}
s1.si = map_surface_index(s1.si);
if(create_it)
{
mesh.AddSegment(s1);
mesh.AddSegment(s2);
}
s1.surfnr1 = s1.surfnr2 = s2.surfnr1 = s2.surfnr2 = -1;
mesh.AddSegment(s1);
mesh.AddSegment(s2);
}
// segi[0] = mapto[segi[0]] not working somehow?
mesh[sei][0] = mapto[segi[0]];
mesh[sei][1] = mapto[segi[1]];
mesh[sej][0] = mapto[segj[0]];
mesh[sej][1] = mapto[segj[1]];
}
// remap the segments to the new points
mesh[sei][0] = mapto[mesh[sei][0]];
mesh[sei][1] = mapto[mesh[sei][1]];
mesh[sej][1] = mapto[mesh[sej][1]];
mesh[sej][0] = mapto[mesh[sej][0]];
}
}
}
else
{
@ -455,9 +447,9 @@ namespace netgen
{
for(SurfaceElementIndex si = 0; si < nse; si++)
{
if(blp.surfid.Contains(mesh[si].GetIndex()))
const auto& sel = mesh[si];
if(blp.surfid.Contains(sel.GetIndex()))
{
const auto& sel = mesh[si];
Element2d newel = sel;
newel.SetIndex(map_surface_index(sel.GetIndex()));
mesh.AddSurfaceElement(newel);
@ -546,36 +538,28 @@ namespace netgen
for(SurfaceElementIndex sei : Range(nse))
{
auto& sel = mesh[sei];
bool to_move = blp.surfid.Contains(sel.GetIndex());
if(blp.domains.Size())
if(!blp.surfid.Contains(sel.GetIndex()))
{
if(blp.outside)
to_move |= blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()];
else
to_move |= !blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()];
}
if(to_move)
{
for(auto& pnum : sel.PNums())
// Check (Doublecheck) if the corresponding point has a
// copy available for remapping
if(mapto[pnum].IsValid())
// Map the surface elements to the new points
pnum = mapto[pnum];
const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
if(blp.outside &&
(!blp.domains[fd.DomainIn()] && !blp.domains[fd.DomainOut()]))
continue;
if(!blp.outside &&
(blp.domains[fd.DomainIn()] || blp.domains[fd.DomainOut()]))
continue;
}
for(auto& pnum : sel.PNums())
if(mapto[pnum].IsValid())
pnum = mapto[pnum];
}
for(ElementIndex ei : Range(ne))
{
auto& el = mesh[ei];
bool to_move = blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()];
if(blp.domains.Size() == 0 || to_move)
// only move the elements on the correct side
if(blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()])
for(auto& pnum : el.PNums())
// Check (Doublecheck) if the corresponding point has a
// copy available for remapping
if(mapto[pnum].IsValid())
// Map the volume elements to the new points
pnum = mapto[pnum];
}

30
libsrc/meshing/curvedelems.cpp
View File

@ -748,21 +748,19 @@ namespace netgen
for (int i2 = 0; i2 < edgenrs.Size(); i2++)
{
// PointIndex pi1 = el[edges[i2][0]];
// PointIndex pi2 = el[edges[i2][1]];
// bool swap = pi1 > pi2;
// Point<3> p1 = mesh[pi1];
// Point<3> p2 = mesh[pi2];
// int order1 = edgeorder[edgenrs[i2]];
// int ndof = max (0, order1-1);
surfnr[edgenrs[i2]] = mesh.GetFaceDescriptor(el.GetIndex()).SurfNr();
gi0[edgenrs[i2]] = el.GeomInfoPi(edges[i2][0]+1);
gi1[edgenrs[i2]] = el.GeomInfoPi(edges[i2][1]+1);
}
auto enr = edgenrs[i2];
surfnr[enr] = mesh.GetFaceDescriptor(el.GetIndex()).SurfNr();
if (el[edges[i2][0]] < el[edges[i2][1]])
{
gi0[enr] = el.GeomInfoPi(edges[i2][0]+1);
gi1[enr] = el.GeomInfoPi(edges[i2][1]+1);
}
else
{
gi1[enr] = el.GeomInfoPi(edges[i2][0]+1);
gi0[enr] = el.GeomInfoPi(edges[i2][1]+1);
}
}
}
@ -1303,6 +1301,8 @@ namespace netgen
SurfaceElementIndex sei = top.GetFace2SurfaceElement (f+1)-1;
if (sei != SurfaceElementIndex(-1)) {
PointGeomInfo gi = mesh[sei].GeomInfoPi(1);
gi.u = 1.0/3.0*(mesh[sei].GeomInfoPi(1).u+mesh[sei].GeomInfoPi(2).u+mesh[sei].GeomInfoPi(3).u);
gi.v = 1.0/3.0*(mesh[sei].GeomInfoPi(1).v+mesh[sei].GeomInfoPi(2).v+mesh[sei].GeomInfoPi(3).v);
geo.ProjectPointGI(surfnr[facenr], pp, gi);
}
else

2
libsrc/meshing/meshclass.cpp
View File

@ -6541,7 +6541,7 @@ namespace netgen
return defaultstring;
if (bcnr < 0 || bcnr >= bcnames.Size())
throw NgException ("illegal bc-number");
throw RangeException("Illegal bc number ", bcnr, 0, bcnames.Size());
if ( bcnames[bcnr] )
return *bcnames[bcnr];

45
libsrc/meshing/meshfunc2d.cpp
View File

@ -20,6 +20,19 @@ namespace netgen
}
mesh.Compress();
bool optimize_swap_separate_faces = false;
if(!mp.quad)
{
bool mixed = false;
ParallelFor( Range(mesh.GetNSE()), [&] (auto i) NETGEN_LAMBDA_INLINE
{
if (mesh[SurfaceElementIndex(i)].GetNP() != 3)
mixed = true;
});
if(mixed)
optimize_swap_separate_faces = true;
}
const char * optstr = mp.optimize2d.c_str();
int optsteps = mp.optsteps2d;
@ -31,16 +44,42 @@ namespace netgen
{
case 's':
{ // topological swap
MeshOptimize2d meshopt(mesh);
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.EdgeSwapping (0);
if(optimize_swap_separate_faces)
{
for(auto i : Range(1, mesh.GetNFD()+1))
{
meshopt.SetFaceIndex(i);
meshopt.EdgeSwapping (0);
}
}
else
{
meshopt.SetFaceIndex(0);
meshopt.EdgeSwapping (0);
}
break;
}
case 'S':
{ // metric swap
MeshOptimize2d meshopt(mesh);
meshopt.SetMetricWeight (mp.elsizeweight);
meshopt.EdgeSwapping (1);
if(optimize_swap_separate_faces)
{
for(auto i : Range(1, mesh.GetNFD()+1))
{
meshopt.SetFaceIndex(i);
meshopt.EdgeSwapping (1);
}
}
else
{
meshopt.SetFaceIndex(0);
meshopt.EdgeSwapping (1);
}
break;
}
case 'm':

2
libsrc/meshing/meshing.hpp
View File

@ -57,10 +57,12 @@ namespace netgen
#include "hprefinement.hpp"
#include "boundarylayer.hpp"
#include "specials.hpp"
}
#include "validate.hpp"
#include "basegeom.hpp"
#include "surfacegeom.hpp"
#ifdef PARALLEL
#include "paralleltop.hpp"

6
libsrc/meshing/parallelmesh.cpp
View File

@ -798,7 +798,7 @@ namespace netgen
// slaves receive the mesh from the master
// workers receive the mesh from the master
void Mesh :: ReceiveParallelMesh ( )
{
int timer = NgProfiler::CreateTimer ("ReceiveParallelMesh");
@ -1058,7 +1058,7 @@ namespace netgen
// distribute the mesh to the slave processors
// distribute the mesh to the worker processors
// call it only for the master !
void Mesh :: Distribute ()
{
@ -1358,7 +1358,7 @@ namespace netgen
// distribute the mesh to the slave processors
// distribute the mesh to the worker processors
// call it only for the master !
void Mesh :: Distribute (NgArray<int> & volume_weights , NgArray<int> & surface_weights, NgArray<int> & segment_weights)
{

57
libsrc/meshing/python_mesh.cpp
View File

@ -949,8 +949,22 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
{
regex pattern(*get_if<string>(&boundary));
for(int i = 1; i<=self.GetNFD(); i++)
if(regex_match(self.GetFaceDescriptor(i).GetBCName(), pattern))
blp.surfid.Append(i);
{
auto& fd = self.GetFaceDescriptor(i);
if(regex_match(fd.GetBCName(), pattern))
{
auto dom_pattern = get_if<string>(&domain);
// only add if adjacent to domain
if(dom_pattern)
{
regex pattern(*dom_pattern);
if(regex_match(self.GetMaterial(fd.DomainIn()), pattern) || (fd.DomainOut() > 0 ? regex_match(self.GetMaterial(fd.DomainOut()), pattern) : false))
blp.surfid.Append(i);
}
else
blp.surfid.Append(i);
}
}
}
if(double* pthickness = get_if<double>(&thickness); pthickness)
@ -1132,6 +1146,45 @@ grow_edges : bool = False
m.def("ReadCGNSFile", &ReadCGNSFile, py::arg("filename"), py::arg("base")=1, "Read mesh and solution vectors from CGNS file");
py::class_<SurfaceGeometry, NetgenGeometry, shared_ptr<SurfaceGeometry>> (m, "SurfaceGeometry")
.def(py::init<>())
.def(py::init([](py::object pyfunc)
{
std::function<Vec<3> (Point<2>)> func = [pyfunc](Point<2> p)
{
py::gil_scoped_acquire aq;
py::tuple pyres = py::extract<py::tuple>(pyfunc(p[0],p[1],0.0)) ();
return Vec<3>(py::extract<double>(pyres[0])(),py::extract<double>(pyres[1])(),py::extract<double>(pyres[2])());
};
auto geo = make_shared<SurfaceGeometry>(func);
return geo;
}), py::arg("mapping"))
.def(NGSPickle<SurfaceGeometry>())
.def("GenerateMesh", [](shared_ptr<SurfaceGeometry> geo,
bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames)
{
if (py::len(py_bbbpts) != py::len(py_bbbnames))
throw Exception("In SurfaceGeometry::GenerateMesh bbbpts and bbbnames do not have same lengths.");
Array<Point<3>> bbbpts(py::len(py_bbbpts));
Array<string> bbbname(py::len(py_bbbpts));
for(int i = 0; i<py::len(py_bbbpts);i++)
{
py::tuple pnt = py::extract<py::tuple>(py_bbbpts[i])();
bbbpts[i] = Point<3>(py::extract<double>(pnt[0])(),py::extract<double>(pnt[1])(),py::extract<double>(pnt[2])());
bbbname[i] = py::extract<string>(py_bbbnames[i])();
}
auto mesh = make_shared<Mesh>();
SetGlobalMesh (mesh);
mesh->SetGeometry(geo);
ng_geometry = geo;
auto result = geo->GenerateMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname);
if(result != 0)
throw Exception("SurfaceGeometry: Meshing failed!");
return mesh;
}, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list())
;
;
}
PYBIND11_MODULE(libmesh, m) {

419
libsrc/meshing/surfacegeom.cpp Normal file
View File

@ -0,0 +1,419 @@
/* *************************************************************************/
/* File: surfacegeom.cpp */
/* Author: Michael Neunteufel */
/* Date: Jun. 2020 */
/* *************************************************************************/
#include <meshing.hpp>
namespace netgen
{
SurfaceGeometry :: SurfaceGeometry()
{
//identity
func = [](Point<2> p) { return Vec<3>(p[0],p[1],0.0); };
}
SurfaceGeometry :: SurfaceGeometry(function<Vec<3>(Point<2>)> _func) : func(_func)
{
;
}
SurfaceGeometry :: SurfaceGeometry(const SurfaceGeometry& geom) : func(geom.func), eps(geom.eps)
{
;
}
void SurfaceGeometry :: CalcHesse(double u, double v, Vec<3>& f_uu, Vec<3>& f_vv, Vec<3>& f_uv) const
{
Point<2> p = Point<2>(u,v);
double pr = p[0]+eps;
double pl = p[0]-eps;
double prr = p[0]+2*eps;
double pll = p[0]-2*eps;
auto dr = GetTangentVectors( pr, v );
auto dl = GetTangentVectors( pl, v );
auto drr = GetTangentVectors( prr, v );
auto dll = GetTangentVectors( pll, v );
f_uu = (1.0/(12.0*eps)) * (8.0*dr[0]-8.0*dl[0]-drr[0]+dll[0]);
f_uv = (1.0/(12.0*eps)) * (8.0*dr[1]-8.0*dl[1]-drr[1]+dll[1]);
pr = p[1]+eps;
pl = p[1]-eps;
prr = p[1]+2*eps;
pll = p[1]-2*eps;
dr = GetTangentVectors(u, pr);
dl = GetTangentVectors(u, pl);
drr = GetTangentVectors(u, prr);
dll = GetTangentVectors(u, pll);
f_vv = (1.0/(12.0*eps)) * (8.0*dr[1]-8.0*dl[1]-drr[1]+dll[1]);
}
Array<Vec<3>> SurfaceGeometry :: GetTangentVectors(double u, double v) const
{
Array<Vec<3>> tang(2);
Point<2> pru = Point<2>(u+eps,v);
Point<2> plu = Point<2>(u-eps,v);
Point<2> prru = Point<2>(u+2*eps,v);
Point<2> pllu = Point<2>(u-2*eps,v);
Point<2> prv = Point<2>(u,v+eps);
Point<2> plv = Point<2>(u,v-eps);
Point<2> prrv = Point<2>(u,v+2*eps);
Point<2> pllv = Point<2>(u,v-2*eps);
tang[0] = 1/(12.0*eps)*( 8.0*func(pru) - 8.0*func(plu) - func(prru) + func(pllu) );
tang[1] = 1/(12.0*eps)*( 8.0*func(prv) - 8.0*func(plv) - func(prrv) + func(pllv) );
return tang;
}
Vec<3> SurfaceGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const
{
Array<Vec<3>> tang = GetTangentVectors(gi->u, gi->v);
auto normal = Cross(tang[0], tang[1]);
return Cross(tang[0], tang[1]);
}
PointGeomInfo SurfaceGeometry :: ProjectPoint(int surfind, Point<3> & p) const
{
throw Exception("In SurfaceGeometry::ProjectPoint");
}
void SurfaceGeometry :: ProjectPointEdge (int surfind, int surfind2, Point<3> & p,
EdgePointGeomInfo* gi) const
{
if (gi == nullptr)
throw Exception("In SurfaceGeometry::ProjectPointEdge: gi is nullptr");
throw Exception("In SurfaceGeometry::ProjectPointEdge: not implemented");
}
bool SurfaceGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
{
Array<Vec<3>> tangs;
Vec<3> diff, f_uu, f_vv, f_uv;
Vec<2> r, dx;
double norm_r, det, energy=0.0, new_energy=0.0, alpha=2.0,u=0.0,v=0.0;
Mat<2,2> mat, inv;
int num=0, maxit=20;
double damping=0.2;
//Solve minimization problem
// argmin_(u,v) 0.5*\| f(u,v)-p\|^2
//via Neton's method:
// F(u,v) = ( (f(u,v)-p)*f_u(u,v), (f(u,v)-p)*f_v(u,v))^T = (0,0)^T
//Stiffness matrix
// F'(u,v) = ( f_u*f_u + (f-p)*f_uu, f_v*f_u + (f-p)*f_uv, f_v*f_u + (f-p)*f_uv, f_v*f_v + (f-p)*f_vv )
do
{
num++;
tangs = GetTangentVectors(gi.u, gi.v);
diff = func(Point<2>(gi.u, gi.v)) - Vec<3>(p);
energy = diff.Length2();
r = Vec<2>( diff*tangs[0], diff*tangs[1] );
norm_r = r.Length2();
CalcHesse(gi.u, gi.v, f_uu, f_vv, f_uv);
mat(0,0) = tangs[0]*tangs[0] + diff*f_uu;
mat(1,0) = mat(0,1) = tangs[0]*tangs[1]+diff*f_uv;
mat(1,1) = tangs[1]*tangs[1]+diff*f_vv;
CalcInverse(mat,inv);
dx = inv*r;
//Linesearch
alpha = 2.0;
do
{
alpha /= 2.0;
u = gi.u - min(1.0,alpha*damping*num)*dx[0];
v = gi.v - min(1.0,alpha*damping*num)*dx[1];
diff = func(Point<2>(u, v)) - Vec<3>(p);
new_energy = diff.Length2();
}
while (alpha > 1e-10 && new_energy > energy+1e-14);
if (alpha <= 1e-10)
throw Exception("In SurfaceGeometry::ProjectPointGI: Linesearch min alpha reached!");
gi.u = u;
gi.v = v;
}
while ( norm_r > 1e-12 && num < maxit);
//Stay in reference domain [0,1]^2
if (gi.u < 0 || gi.u > 1 || gi.v < 0 || gi.v > 1)
{
cout << "Warning: Projected point outside [0,1]^2: u=" << gi.u << ",v=" << gi.v <<". Setting back." << endl;
gi.u = min(max(gi.u,0.0),1.0);
gi.v = min(max(gi.v,0.0),1.0);
}
p = Point<3>(func(Point<2>(gi.u,gi.v)));
if (num == maxit)
{
//cout << "In SurfaceGeometry::ProjectPointGI: Newton did not converge" << endl;
throw Exception("In SurfaceGeometry::ProjectPointGI: Newton did not converge");
}
return true;
}
bool SurfaceGeometry :: CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const
{
throw Exception("In SurfaceGeometry::CalcPointGeomInfo: not implemented");
return false;
}
void SurfaceGeometry :: 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
{
newp = p1+secpoint*(p2-p1);
PointGeomInfo pgi;
pgi.u = ap1.u+secpoint*(ap2.u-ap1.u);
pgi.v = ap1.v+secpoint*(ap2.v-ap1.v);
ProjectPointGI(surfi1, newp, pgi);
newgi.u = pgi.u;
newgi.v = pgi.v;
newgi.edgenr = ap1.edgenr;
newgi.body = -1;
newgi.dist = -1.0;
}
void SurfaceGeometry :: 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.u = gi1.u+secpoint*(gi2.u-gi1.u);
newgi.v = gi1.v+secpoint*(gi2.v-gi1.v);
newgi.trignum = -1;
ProjectPointGI(surfi, newp, newgi);
}
int SurfaceGeometry :: GenerateMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames)
{
mesh->SetDimension(3);
Array<bool> found(bbbpts.Size());
found = false;
Array<PointIndex> indbbbpts(bbbpts.Size());
Array<PointIndex> pids;
Array<PointGeomInfo> pgis;
for(int i=0; i <= ny; i++)
for(int j=0; j <= nx; j++)
{
PointGeomInfo pgi;
pgi.trignum = -1;
pgi.u = double(j)/nx;
pgi.v = double(i)/ny;
Point<3> pnt = Point<3>(func(Point<2>(pgi.u,pgi.v)));
pids.Append(mesh->AddPoint(pnt));
pgis.Append(pgi);
for (int k = 0; k < bbbpts.Size(); k++)
{
auto diff = pnt - bbbpts[k];
if(diff.Length2() < 1e-14)
{
found[k] = true;
indbbbpts[k] = pids[pids.Size()-1];
}
}
}
for (bool f : found)
if (!f)
throw Exception("In SurfaceGeometry :: GenerateMesh: bbbpts not resolved in mesh.");
FaceDescriptor fd;
fd.SetSurfNr(1);
fd.SetDomainIn(1);
fd.SetDomainOut(0);
fd.SetBCProperty(1);
mesh->AddFaceDescriptor(fd);
for(int i=0; i < ny; i++)
{
for(int j=0; j < nx; j++)
{
int base = i * (nx+1) + j;
if (quads)
{
int pnum[4] = {base,base+1,base+nx+2,base+nx+1};
Element2d el = Element2d(QUAD);
for (int i = 0; i < 4; i++)
{
el[i] = pids[pnum[i]];
el.GeomInfoPi(i+1) = pgis[pnum[i]];
}
el.SetIndex(1);
mesh->AddSurfaceElement(el);
}
else
{
Array<int> pnum1(3);
Array<int> pnum2(3);
if (flip_triangles)
{
pnum1[0] = base;
pnum1[1] = base+1;
pnum1[2] = base+nx+2;
pnum2[0] = base;
pnum2[1] = base+nx+2;
pnum2[2] = base+nx+1;
}
else
{
pnum1[0] = base;
pnum1[1] = base+1;
pnum1[2] = base+nx+1;
pnum2[0] = base+1;
pnum2[1] = base+nx+2;
pnum2[2] = base+nx+1;
}
Element2d el = Element2d(TRIG);
for (int i = 0; i < 3; i++)
{
el[i] = pids[pnum1[i]];
el.GeomInfoPi(i+1) = pgis[pnum1[i]];
}
el.SetIndex(1);
mesh->AddSurfaceElement(el);
for (int i = 0; i < 3; i++)
{
el[i] = pids[pnum2[i]];
el.GeomInfoPi(i+1) = pgis[pnum2[i]];
}
mesh->AddSurfaceElement(el);
}
}
}
Segment seg;
seg.si = 1;
seg.edgenr = 1;
seg.epgeominfo[0].edgenr = 1;
seg.epgeominfo[1].edgenr = 1;
// needed for codim2 in 3d
seg.edgenr = 1;
for(int i=0; i < nx; i++)
{
seg[0] = pids[i];
seg[1] = pids[i+1];
seg.geominfo[0] = pgis[i];
seg.geominfo[1] = pgis[i+1];
seg.epgeominfo[0].u = pgis[i].u;
seg.epgeominfo[0].v = pgis[i].v;
seg.epgeominfo[0].edgenr = seg.edgenr;
seg.epgeominfo[1].u = pgis[i+1].u;
seg.epgeominfo[1].v = pgis[i+1].v;
seg.epgeominfo[1].edgenr = seg.edgenr;
mesh->AddSegment(seg);
}
seg.si = 2;
seg.edgenr = 2;
for(int i=0; i<ny; i++)
{
seg[0] = pids[i*(nx+1)+nx];
seg[1] = pids[(i+1)*(nx+1)+nx];
seg.geominfo[0] = pgis[i*(nx+1)+nx];
seg.geominfo[1] = pgis[(i+1)*(nx+1)+nx];
seg.epgeominfo[0].u = pgis[i*(nx+1)+nx].u;
seg.epgeominfo[0].v = pgis[i*(nx+1)+nx].v;
seg.epgeominfo[0].edgenr = seg.edgenr;
seg.epgeominfo[1].u = pgis[(i+1)*(nx+1)+nx].u;
seg.epgeominfo[1].v = pgis[(i+1)*(nx+1)+nx].v;
seg.epgeominfo[1].edgenr = seg.edgenr;
mesh->AddSegment(seg);
}
seg.si = 3;
seg.edgenr = 3;
for(int i=0; i<nx; i++)
{
seg[0] = pids[ny*(nx+1)+i+1];
seg[1] = pids[ny*(nx+1)+i];
seg.geominfo[0] = pgis[ny*(nx+1)+i+1];
seg.geominfo[1] = pgis[ny*(nx+1)+i];
seg.epgeominfo[0].u = pgis[ny*(nx+1)+i+1].u;
seg.epgeominfo[0].v = pgis[ny*(nx+1)+i+1].v;
seg.epgeominfo[0].edgenr = seg.edgenr;
seg.epgeominfo[1].u = pgis[ny*(nx+1)+i].u;
seg.epgeominfo[1].v = pgis[ny*(nx+1)+i].v;
seg.epgeominfo[1].edgenr = seg.edgenr;
mesh->AddSegment(seg);
}
seg.si = 4;
seg.edgenr = 4;
for(int i=0; i<ny; i++)
{
seg[0] = pids[(i+1)*(nx+1)];
seg[1] = pids[i*(nx+1)];
seg.geominfo[0] = pgis[(i+1)*(nx+1)];
seg.geominfo[1] = pgis[i*(nx+1)];
seg.epgeominfo[0].u = pgis[(i+1)*(nx+1)].u;
seg.epgeominfo[0].v = pgis[(i+1)*(nx+1)].v;
seg.epgeominfo[0].edgenr = seg.edgenr;
seg.epgeominfo[1].u = pgis[i*(nx+1)].u;
seg.epgeominfo[1].v = pgis[i*(nx+1)].v;
seg.epgeominfo[1].edgenr = seg.edgenr;
mesh->AddSegment(seg);
}
mesh->SetCD2Name(1, "bottom");
mesh->SetCD2Name(2, "right");
mesh->SetCD2Name(3, "top");
mesh->SetCD2Name(4, "left");
for (int i = 0; i < bbbpts.Size(); i++)
{
Element0d el;
el.pnum = indbbbpts[i];
el.index = i+1;
mesh->pointelements.Append(el);
mesh->SetCD3Name(i+1, bbbnames[i]);
}
mesh->Compress();
mesh->UpdateTopology();
return 0;
}
};

70
libsrc/meshing/surfacegeom.hpp Normal file
View File

@ -0,0 +1,70 @@
#ifndef FILE_SURFACEGEOM
#define FILE_SURFACEGEOM
/* *************************************************************************/
/* File: surfacegeom.hpp */
/* Author: Michael Neunteufel */
/* Date: Jun. 2020 */
/* *************************************************************************/
#include <functional>
namespace netgen
{
class DLL_HEADER SurfaceGeometry : public NetgenGeometry
{
function<Vec<3>(Point<2>)> func;
double eps=1e-4;
private:
void CalcHesse(double u, double v, Vec<3>& f_uu, Vec<3>& f_vv, Vec<3>& f_uv) const;
public:
SurfaceGeometry();
SurfaceGeometry(function<Vec<3>(Point<2>)> func);
SurfaceGeometry(const SurfaceGeometry& geom);
SurfaceGeometry& operator =(const SurfaceGeometry& geom)
{
func = geom.func;
eps = geom.eps;
return *this;
}
Array<Vec<3>> GetTangentVectors(double u, double v) const;
virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override;
virtual PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override;
virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p,
EdgePointGeomInfo* gi = nullptr) const override;
virtual bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
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 override;
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;
int GenerateMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames);
};
}
#endif //SURFACEGEOM

83
libsrc/occ/occgeom.cpp
View File

@ -169,8 +169,91 @@ namespace netgen
}
void OCCGeometry :: GlueGeometry()
{
PrintMessage(1, "OCC Glue Geometry");
/*
//
BRep_Builder builder;
TopoDS_Shape my_fuse;
int cnt = 0;
for (TopExp_Explorer exp_solid(shape, TopAbs_SOLID); exp_solid.More(); exp_solid.Next())
{
cout << "cnt = " << cnt << endl;
if (cnt == 0)
my_fuse = exp_solid.Current();
else
// my_fuse = BRepAlgoAPI_Fuse (my_fuse, exp_solid.Current());
my_fuse = QANewModTopOpe_Glue::QANewModTopOpe_Glue(my_fuse, exp_solid.Current());
cnt++;
}
cout << "remove" << endl;
// for (int i = 1; i <= somap.Size(); i++)
// builder.Remove (shape, somap(i));
cout << "now add" << endl;
// builder.Add (shape, my_fuse);
shape = my_fuse;
cout << "build fmap" << endl;
BuildFMap();
*/
// from
// https://www.opencascade.com/doc/occt-7.4.0/overview/html/occt_user_guides__boolean_operations.html
BOPAlgo_Builder aBuilder;
// Setting arguments
TopTools_ListOfShape aLSObjects;
for (TopExp_Explorer exp_solid(shape, TopAbs_SOLID); exp_solid.More(); exp_solid.Next())
aLSObjects.Append (exp_solid.Current());
aBuilder.SetArguments(aLSObjects);
// Setting options for GF
// Set parallel processing mode (default is false)
// Standard_Boolean bRunParallel = Standard_True;
// aBuilder.SetRunParallel(bRunParallel);
// Set Fuzzy value (default is Precision::Confusion())
// Standard_Real aFuzzyValue = 1.e-5;
// aBuilder.SetFuzzyValue(aFuzzyValue);
// Set safe processing mode (default is false)
// Standard_Boolean bSafeMode = Standard_True;
// aBuilder.SetNonDestructive(bSafeMode);
// Set Gluing mode for coinciding arguments (default is off)
// BOPAlgo_GlueEnum aGlue = BOPAlgo_GlueShift;
// aBuilder.SetGlue(aGlue);
// Disabling/Enabling the check for inverted solids (default is true)
// Standard Boolean bCheckInverted = Standard_False;
// aBuilder.SetCheckInverted(bCheckInverted);
// Set OBB usage (default is false)
// Standard_Boolean bUseOBB = Standard_True;
// aBuilder.SetUseOBB(buseobb);
// Perform the operation
aBuilder.Perform();
// Check for the errors
#if OCC_VERSION_HEX >= 0x070200
if (aBuilder.HasErrors())
{
cout << "builder has errors" << endl;
return;
}
// Check for the warnings
if (aBuilder.HasWarnings())
{
// treatment of the warnings
;
}
#endif
// result of the operation
shape = aBuilder.Shape();
BuildFMap();
}
void OCCGeometry :: HealGeometry ()
{
int nrc = 0, nrcs = 0,

3
libsrc/occ/occgeom.hpp
View File

@ -78,7 +78,7 @@
#include "Bnd_Box.hxx"
#include "ShapeAnalysis.hxx"
#include "ShapeBuild_ReShape.hxx"
#include "BOPAlgo_Builder.hxx"
// Philippose - 29/01/2009
// OpenCascade XDE Support
@ -343,6 +343,7 @@ namespace netgen
void MakeSolid();
void HealGeometry();
void GlueGeometry();
// Philippose - 15/01/2009
// Sets the maximum mesh size for a given face

1
libsrc/occ/python_occ.cpp
View File

@ -66,6 +66,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
}), py::arg("filename"),
"Load OCC geometry from step, brep or iges file")
.def(NGSPickle<OCCGeometry>())
.def("Glue", &OCCGeometry::GlueGeometry)
.def("Heal",[](OCCGeometry & self, double tolerance, bool fixsmalledges, bool fixspotstripfaces, bool sewfaces, bool makesolids, bool splitpartitions)
{
self.tolerance = tolerance;

22
tests/pytest/test_splinegeo_tensordomainmeshing.py Normal file
View File

@ -0,0 +1,22 @@
from netgen.geom2d import *
def test_tensordomainmeshing():
geo = SplineGeometry()
w = 10
h = 0.01
p = [ (0, 0), (w, 0), (w, h), (0, h) ]
p = [geo.AppendPoint(*px) for px in p]
l0 = geo.Append ( ["line", p[0], p[1]], leftdomain=1, rightdomain=0 )
l1 = geo.Append ( ["line", p[1], p[2]], leftdomain=1, rightdomain=0)
geo.Append ( ["line", p[3], p[2]], leftdomain=0, rightdomain=1, copy=l0 )
geo.Append ( ["line", p[0], p[3]], leftdomain=0, rightdomain=1, copy=l1 )
geo._SetDomainTensorMeshing(1, True)
mesh = geo.GenerateMesh(maxh=1)
for el in mesh.Elements2D():
print(el.vertices)
assert len(el.vertices) == 4