ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-26 05:50:32 +05:00
Code Issues Packages Projects Releases Wiki Activity

merge master

Browse Source
This commit is contained in:
Christopher Lackner 2016-10-31 15:37:02 +01:00
commit 2d77b35b39
13 changed files with 191 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
CMakeLists.txt
View File

@ -364,11 +364,7 @@ if(UNIX)
set(CPACK_PACKAGE_NAME "${CPACK_PACKAGE_NAME}_mpi")
endif(USE_MPI)
if(USE_OCC)
if("${UBUNTU_VERSION}" STREQUAL "xenial")
set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf10")
else()
set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf8")
endif()
set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf-dev")
endif(USE_OCC)
set(CPACK_DEBIAN_PACKAGE_SECTION Science)
set(CPACK_DEBIAN_PACKAGE_NAME ${CPACK_PACKAGE_NAME})

19
libsrc/csg/python_csg.cpp
View File

@ -167,6 +167,9 @@ namespace netgen
extern CSGeometry * ParseCSG (istream & istr);
}
static Transformation<3> global_trafo(Vec<3> (0,0,0));
DLL_HEADER void ExportCSG()
{
ModuleScope module("csg");
@ -186,10 +189,22 @@ DLL_HEADER void ExportCSG()
;
bp::def ("Pnt", FunctionPointer
([](double x, double y, double z) { return Point<3>(x,y,z); }));
([](double x, double y, double z) { return global_trafo(Point<3>(x,y,z)); }));
bp::def ("Pnt", FunctionPointer
([](double x, double y) { return Point<2>(x,y); }));
bp::def ("SetTransformation", FunctionPointer
([](int dir, double angle)
{
if (dir > 0)
global_trafo.SetAxisRotation (dir, angle*M_PI/180);
else
global_trafo = Transformation<3> (Vec<3>(0,0,0));
}),
(bp::arg("dir")=int(0), bp::arg("angle")=0));
bp::class_<Vec<2>> ("Vec2d", bp::init<double,double>())
.def ("__str__", &ToString<Vec<3>>)
.def(bp::self+bp::self)
@ -209,7 +224,7 @@ DLL_HEADER void ExportCSG()
;
bp::def ("Vec", FunctionPointer
([] (double x, double y, double z) { return Vec<3>(x,y,z); }));
([] (double x, double y, double z) { return global_trafo(Vec<3>(x,y,z)); }));
bp::def ("Vec", FunctionPointer
([] (double x, double y) { return Vec<2>(x,y); }));

2
libsrc/general/hashtabl.hpp
View File

@ -779,7 +779,7 @@ protected:
mask = size-1;
// cout << "mask = " << mask << endl;
invalid = -1;
for (int i = 0; i < size; i++)
for (size_t i = 0; i < size; i++)
hash[i].I1() = invalid;
}

16
libsrc/general/mysimd.hpp
View File

@ -54,8 +54,18 @@ namespace netgen
#ifdef __AVX__
template <typename T>
class AlignedAlloc
{
public:
void * operator new (size_t s) { return _mm_malloc(s, alignof(T)); }
void * operator new[] (size_t s) { return _mm_malloc(s, alignof(T)); }
void operator delete (void * p) { _mm_free(p); }
void operator delete[] (void * p) { _mm_free(p); }
};
template<>
class alignas(32) SIMD<double>
class alignas(32) SIMD<double> : public AlignedAlloc<SIMD<double>>
{
__m256d data;
@ -158,6 +168,10 @@ namespace netgen
#else
// it's only a dummy without AVX
template <typename T>
class AlignedAlloc { ; };
template<>
class SIMD<double>
{

9
libsrc/geom2d/python_geom2d.cpp
View File

@ -51,7 +51,8 @@ DLL_HEADER void ExportGeom2d()
return self.geompoints.Size()-1;
}),
(bp::arg("self"), bp::arg("x"), bp::arg("y"), bp::arg("maxh") = 1e99, bp::arg("hpref")=false))
.def("Append", FunctionPointer([](SplineGeometry2d &self, bp::list segment, int leftdomain, int rightdomain, bp::object bc, double maxh, bool hpref)
.def("Append", FunctionPointer([](SplineGeometry2d &self, bp::list segment, int leftdomain, int rightdomain,
bp::object bc, bp::object copy, double maxh, bool hpref)
{
bp::extract<std::string> segtype(segment[0]);
@ -85,6 +86,9 @@ DLL_HEADER void ExportGeom2d()
seg->hpref_right = hpref;
seg->reffak = 1;
seg->copyfrom = -1;
if (bp::extract<int>(copy).check())
seg->copyfrom = bp::extract<int>(copy)()+1;
if (bp::extract<int>(bc).check())
seg->bc = bp::extract<int>(bc)();
else if (bp::extract<string>(bc).check())
@ -98,8 +102,9 @@ DLL_HEADER void ExportGeom2d()
else
seg->bc = self.GetNSplines()+1;
self.AppendSegment(seg);
return self.GetNSplines()-1;
}), (bp::arg("self"), bp::arg("point_indices"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0,
bp::arg("bc")=bp::object(), bp::arg("maxh")=1e99, bp::arg("hpref")=false
bp::arg("bc")=bp::object(), bp::arg("copy")=bp::object(), bp::arg("maxh")=1e99, bp::arg("hpref")=false
))

5
libsrc/gprim/geomfuncs.hpp
View File

@ -82,9 +82,10 @@ namespace netgen
*/
// inline Vec<3> Cross (const Vec<3> & v1, const Vec<3> & v2)
inline Vec<3> Cross (Vec<3> v1, Vec<3> v2)
template <typename T>
inline Vec<3,T> Cross (Vec<3,T> v1, Vec<3,T> v2)
{
return Vec<3>
return Vec<3,T>
( v1(1) * v2(2) - v1(2) * v2(1),
v1(2) * v2(0) - v1(0) * v2(2),
v1(0) * v2(1) - v1(1) * v2(0) );

25
libsrc/gprim/geomobjects.hpp
View File

@ -17,7 +17,7 @@ namespace netgen
template <int D, typename T>
class Point
class Point : public AlignedAlloc<Point<D,T>>
{
protected:
@ -39,8 +39,9 @@ namespace netgen
Point (T ax, T ay, T az, T au)
{ x[0] = ax; x[1] = ay; x[2] = az; x[3] = au;}
Point (const Point<D> & p2)
{ for (int i = 0; i < D; i++) x[i] = p2.x[i]; }
template <typename T2>
Point (const Point<D,T2> & p2)
{ for (int i = 0; i < D; i++) x[i] = p2(i); }
explicit Point (const Vec<D> & v)
{ for (int i = 0; i < D; i++) x[i] = v(i); }
@ -65,7 +66,7 @@ namespace netgen
};
template <int D, typename T>
class Vec
class Vec : public AlignedAlloc<Vec<D,T>>
{
protected:
@ -90,15 +91,15 @@ namespace netgen
Vec (const Vec<D> & p2)
{ for (int i = 0; i < D; i++) x[i] = p2.x[i]; }
explicit Vec (const Point<D> & p)
explicit Vec (const Point<D,T> & p)
{ for (int i = 0; i < D; i++) x[i] = p(i); }
Vec (const Vec<D> & p1, const Vec<D> & p2)
Vec (const Vec & p1, const Vec & p2)
{ for(int i=0; i<D; i++) x[i] = p2(i)-p1(1); }
Vec & operator= (const Vec<D> & p2)
Vec & operator= (const Vec & p2)
{
for (int i = 0; i < D; i++) x[i] = p2.x[i];
return *this;
@ -131,12 +132,12 @@ namespace netgen
return l;
}
const Vec<D> & Normalize ()
Vec & Normalize ()
{
T l = Length();
if (l != 0)
for (int i = 0; i < D; i++)
x[i] /= l;
// if (l != 0)
for (int i = 0; i < D; i++)
x[i] /= (l+1e-40);
return *this;
}
@ -148,7 +149,7 @@ namespace netgen
template <int H, int W=H, typename T = double>
class Mat
class Mat : public AlignedAlloc<Mat<H,W,T>>
{
protected:

18
libsrc/gprim/geomops.hpp
View File

@ -29,10 +29,10 @@ namespace netgen
template <int D>
inline Point<D> operator+ (const Point<D> & a, const Vec<D> & b)
template <int D, typename T>
inline Point<D,T> operator+ (const Point<D,T> & a, const Vec<D,T> & b)
{
Point<D> res;
Point<D,T> res;
for (int i = 0; i < D; i++)
res(i) = a(i) + b(i);
return res;
@ -40,10 +40,10 @@ namespace netgen
template <int D>
inline Vec<D> operator- (const Point<D> & a, const Point<D> & b)
template <int D, typename T>
inline Vec<D,T> operator- (const Point<D,T> & a, const Point<D,T> & b)
{
Vec<D> res;
Vec<D,T> res;
for (int i = 0; i < D; i++)
res(i) = a(i) - b(i);
return res;
@ -69,10 +69,10 @@ namespace netgen
template <int D>
inline Vec<D> operator* (double s, const Vec<D> & b)
template <int D, typename T>
inline Vec<D,T> operator* (T s, const Vec<D,T> & b)
{
Vec<D> res;
Vec<D,T> res;
for (int i = 0; i < D; i++)
res(i) = s * b(i);
return res;

3
libsrc/gprim/transform3d.hpp
View File

@ -182,6 +182,9 @@ public:
{
to = m * from;
}
Point<D> operator() (Point<D> from) const { Point<D> to; Transform(from, to); return to; }
Vec<D> operator() (Vec<D> from) const { Vec<D> to; Transform(from, to); return to; }
};
template <int D>

9
libsrc/meshing/meshtype.hpp
View File

@ -1024,13 +1024,6 @@ namespace netgen
class DLL_HEADER MeshingParameters
{
public:
/**
3d optimization strategy:
// m .. move nodes
@ -1135,6 +1128,8 @@ namespace netgen
///
MeshingParameters (const MeshingParameters & mp2) = default;
MeshingParameters (MeshingParameters && mp2) = default;
MeshingParameters & operator= (const MeshingParameters & mp2) = default;
MeshingParameters & operator= (MeshingParameters && mp2) = default;
///
void Print (ostream & ost) const;
///

13
libsrc/meshing/python_mesh.cpp
View File

@ -477,14 +477,21 @@ DLL_HEADER void ExportNetgenMeshing()
.def ("GetNCD2Names", &Mesh::GetNCD2Names)
.def ("GenerateVolumeMesh", FunctionPointer
([](Mesh & self)
([](Mesh & self, bp::object pymp)
{
cout << "generate vol mesh" << endl;
MeshingParameters mp;
mp.optsteps3d = 5;
if (bp::extract<MeshingParameters>(pymp).check())
mp = bp::extract<MeshingParameters>(pymp)();
else
{
mp.optsteps3d = 5;
}
MeshVolume (mp, self);
OptimizeVolume (mp, self);
}))
}),
(bp::arg("self"), bp::arg("mp")=bp::object())
)
.def ("OptimizeVolumeMesh", FunctionPointer
([](Mesh & self)

12
libsrc/visualization/soldata.hpp
View File

@ -92,6 +92,18 @@ namespace netgen
return res;
}
#ifdef __AVX__
virtual bool GetMultiSurfValue (size_t selnr, size_t facetnr, size_t npts,
const __m256d * xref,
const __m256d * x,
const __m256d * dxdxref,
__m256d * values)
{
cerr << "GetMultiSurfVaue not overloaded" << endl;
return false;
}
#endif
virtual bool GetSegmentValue (int segnr, double xref, double * values)
{ return false; }

96
libsrc/visualization/vssolution.cpp
View File

@ -1256,6 +1256,14 @@ namespace netgen
Array<double> values(npt);
Array<double> mvalues(npt);
#ifdef __AVX__
Array<Point<2,SIMD<double>> > simd_pref ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
Array<Point<3,SIMD<double>> > simd_points ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
Array<Mat<3,2,SIMD<double>> > simd_dxdxis ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
Array<Vec<3,SIMD<double>> > simd_nvs( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
Array<SIMD<double>> simd_values( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() * sol->components);
#endif
if (sol && sol->draw_surface) mvalues.SetSize (npt * sol->components);
Array<complex<double> > valuesc(npt);
@ -1436,6 +1444,93 @@ namespace netgen
if ( el.GetType() == TRIG || el.GetType() == TRIG6 )
{
#ifdef __AVX_try_it_out__
bool curved = curv.IsSurfaceElementCurved(sei);
for (int iy = 0, ii = 0; iy <= n; iy++)
for (int ix = 0; ix <= n-iy; ix++, ii++)
pref[ii] = Point<2> (ix*invn, iy*invn);
constexpr size_t simd_size = SIMD<double>::Size();
size_t simd_npt = (npt+simd_size-1)/simd_size;
for (size_t i = 0; i < simd_npt; i++)
{
simd_pref[i](0).SIMD_function ([&] (size_t j) { size_t ii = i*simd_size+j; return (ii < npt) ? pref[ii](0) : 0; }, std::true_type());
simd_pref[i](1).SIMD_function ([&] (size_t j) { size_t ii = i*simd_size+j; return (ii < npt) ? pref[ii](1) : 0; }, std::true_type());
}
if (curved)
{
mesh->GetCurvedElements().
CalcMultiPointSurfaceTransformation<3> (sei, simd_npt,
&simd_pref[0](0), 2,
&simd_points[0](0), 3,
&simd_dxdxis[0](0,0), 6);
for (size_t ii = 0; ii < simd_npt; ii++)
simd_nvs[ii] = Cross (simd_dxdxis[ii].Col(0), simd_dxdxis[ii].Col(1)).Normalize();
}
else
{
Point<3,SIMD<double>> p1 = mesh->Point (el[0]);
Point<3,SIMD<double>> p2 = mesh->Point (el[1]);
Point<3,SIMD<double>> p3 = mesh->Point (el[2]);
Vec<3,SIMD<double>> vx = p1-p3;
Vec<3,SIMD<double>> vy = p2-p3;
for (size_t ii = 0; ii < simd_npt; ii++)
{
simd_points[ii] = p3 + simd_pref[ii](0) * vx + simd_pref[ii](1) * vy;
for (size_t j = 0; j < 3; j++)
{
simd_dxdxis[ii](j,0) = vx(j);
simd_dxdxis[ii](j,1) = vy(j);
}
}
Vec<3,SIMD<double>> nv = Cross (vx, vy).Normalize();
for (size_t ii = 0; ii < simd_npt; ii++)
simd_nvs[ii] = nv;
}
bool drawelem = false;
if (sol && sol->draw_surface)
{
drawelem = sol->solclass->GetMultiSurfValue (sei, -1, simd_npt,
&simd_pref[0](0).Data(),
&simd_points[0](0).Data(),
&simd_dxdxis[0](0).Data(),
&simd_values[0].Data());
for (size_t j = 0; j < sol->components; j++)
for (size_t i = 0; i < npt; i++)
mvalues[i*sol->components+j] = ((double*)&simd_values[j*simd_npt])[i];
if (usetexture == 2)
for (int ii = 0; ii < npt; ii++)
valuesc[ii] = ExtractValueComplex(sol, scalcomp, &mvalues[ii*sol->components]);
else
for (int ii = 0; ii < npt; ii++)
values[ii] = ExtractValue(sol, scalcomp, &mvalues[ii*sol->components]);
}
for (size_t i = 0; i < npt; i++)
{
size_t ii = i/4;
size_t r = i%4;
for (int j = 0; j < 2; j++)
pref[i](j) = simd_pref[ii](j)[r];
for (int j = 0; j < 3; j++)
points[i](j) = simd_points[ii](j)[r];
for (int j = 0; j < 3; j++)
nvs[i](j) = simd_nvs[ii](j)[r];
}
if (deform)
for (int ii = 0; ii < npt; ii++)
points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1));
#else
bool curved = curv.IsSurfaceElementCurved(sei);
for (int iy = 0, ii = 0; iy <= n; iy++)
@ -1492,6 +1587,7 @@ namespace netgen
if (deform)
for (int ii = 0; ii < npt; ii++)
points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1));
#endif
int save_usetexture = usetexture;
if (!drawelem)