netgen/libsrc/interface/nginterface_v2.cpp

1112 lines
28 KiB
C++
Raw Normal View History

2011-03-09 20:47:04 +05:00
#include <meshing.hpp>
#ifdef SOCKETS
#include "../sockets/sockets.hpp"
#endif
#include "nginterface.h"
#include "nginterface_v2.hpp"
2015-01-26 19:17:59 +05:00
// #include <visual.hpp>
2011-03-09 20:47:04 +05:00
2016-07-10 21:07:36 +05:00
#include "writeuser.hpp"
2011-03-09 20:47:04 +05:00
namespace netgen
{
2014-09-08 21:21:09 +06:00
extern shared_ptr<Mesh> mesh;
2011-03-09 20:47:04 +05:00
}
namespace netgen
{
2014-01-07 16:42:39 +06:00
#define NGX_INLINE
#include "nginterface_v2_impl.hpp"
2011-03-09 20:47:04 +05:00
2014-09-08 21:21:09 +06:00
shared_ptr<Mesh> Ngx_Mesh :: SelectMesh () const
{
shared_ptr<Mesh> hmesh = netgen::mesh;
2014-10-19 19:53:57 +06:00
2014-09-08 21:21:09 +06:00
netgen::mesh = mesh;
2015-01-09 02:18:33 +05:00
SetGlobalMesh (mesh);
2014-10-19 19:53:57 +06:00
2014-09-08 21:21:09 +06:00
return hmesh;
}
Ngx_Mesh :: Ngx_Mesh (shared_ptr<Mesh> amesh)
2014-01-20 15:30:17 +06:00
{
if (amesh)
mesh = amesh;
else
2014-09-08 21:21:09 +06:00
mesh = netgen::mesh;
2014-01-20 15:30:17 +06:00
}
2014-09-08 21:21:09 +06:00
2013-04-03 02:27:35 +06:00
Ngx_Mesh * LoadMesh (const string & filename)
2012-06-25 22:20:01 +06:00
{
2014-09-08 21:21:09 +06:00
netgen::mesh.reset();
2013-04-03 02:27:35 +06:00
Ng_LoadMesh (filename.c_str());
2014-09-08 21:21:09 +06:00
return new Ngx_Mesh (netgen::mesh);
2012-06-25 22:20:01 +06:00
}
2011-03-09 20:47:04 +05:00
2014-01-07 16:42:39 +06:00
void Ngx_Mesh :: LoadMesh (const string & filename)
{
2014-09-08 21:21:09 +06:00
netgen::mesh.reset();
2014-01-07 16:42:39 +06:00
Ng_LoadMesh (filename.c_str());
2015-07-23 17:11:51 +05:00
// mesh = move(netgen::mesh);
2014-09-08 21:21:09 +06:00
mesh = netgen::mesh;
2014-01-07 16:42:39 +06:00
}
2014-02-24 01:31:40 +06:00
void Ngx_Mesh :: LoadMesh (istream & ist)
{
2014-09-08 21:21:09 +06:00
netgen::mesh = make_shared<Mesh>();
2014-02-24 01:31:40 +06:00
netgen::mesh -> Load (ist);
2015-07-23 17:11:51 +05:00
// mesh = move(netgen::mesh);
2014-09-08 21:21:09 +06:00
mesh = netgen::mesh;
2015-01-09 02:18:33 +05:00
SetGlobalMesh (mesh);
2014-02-24 01:31:40 +06:00
}
void Ngx_Mesh :: SaveMesh (ostream & ost) const
{
mesh -> Save (ost);
}
void Ngx_Mesh :: DoArchive (ngstd::Archive & archive)
{
2018-04-27 11:36:06 +05:00
if (archive.Input()) mesh = make_shared<Mesh>();
mesh->DoArchive(archive);
if (archive.Input())
{
netgen::mesh = mesh;
SetGlobalMesh (mesh);
}
/*
2014-02-24 01:31:40 +06:00
if (archive.Output())
{
stringstream str;
SaveMesh (str);
string st = str.str();
archive & st;
}
else
{
string st;
archive & st;
stringstream str(st);
LoadMesh (str);
}
2018-04-27 11:36:06 +05:00
*/
2014-02-24 01:31:40 +06:00
}
2014-10-19 19:53:57 +06:00
void Ngx_Mesh :: UpdateTopology ()
{
if (mesh)
mesh -> UpdateTopology();
}
2014-02-24 01:31:40 +06:00
2014-01-07 16:42:39 +06:00
/*
2013-04-03 02:27:35 +06:00
Ngx_Mesh :: Ngx_Mesh (Mesh * amesh)
: mesh(amesh)
{ ; }
2014-01-07 16:42:39 +06:00
*/
2013-04-03 02:27:35 +06:00
Ngx_Mesh :: ~Ngx_Mesh ()
2017-08-27 17:51:35 +05:00
{
2017-08-28 20:23:31 +05:00
// causes crashes when global variable netgen::mesh is destructed
// before visualization data
if (mesh == netgen::mesh)
netgen::mesh = nullptr;
2011-03-09 20:47:04 +05:00
}
2013-04-03 02:27:35 +06:00
int Ngx_Mesh :: GetDimension() const
2011-03-09 20:47:04 +05:00
{
2013-04-03 02:27:35 +06:00
return mesh -> GetDimension();
2011-03-09 20:47:04 +05:00
}
2013-04-03 02:27:35 +06:00
int Ngx_Mesh :: GetNLevels() const
2011-03-09 20:47:04 +05:00
{
2013-04-03 02:27:35 +06:00
return mesh -> mglevels;
}
int Ngx_Mesh :: GetNElements (int dim) const
{
switch (dim)
{
2015-11-11 22:46:18 +05:00
case 0: return mesh -> pointelements.Size();
2013-04-03 02:27:35 +06:00
case 1: return mesh -> GetNSeg();
case 2: return mesh -> GetNSE();
case 3: return mesh -> GetNE();
}
return -1;
}
int Ngx_Mesh :: GetNNodes (int nt) const
{
switch (nt)
{
case 0: return mesh -> GetNV();
case 1: return mesh->GetTopology().GetNEdges();
case 2: return mesh->GetTopology().GetNFaces();
case 3: return mesh -> GetNE();
}
return -1;
2011-03-09 20:47:04 +05:00
}
2014-01-07 16:42:39 +06:00
/*
2013-04-03 02:27:35 +06:00
Ng_Point Ngx_Mesh :: GetPoint (int nr) const
{
2014-01-07 16:42:39 +06:00
return Ng_Point (&mesh->Point(nr + PointIndex::BASE)(0));
2013-04-03 02:27:35 +06:00
}
2014-01-07 16:42:39 +06:00
*/
2011-03-09 20:47:04 +05:00
2015-11-11 22:46:18 +05:00
/*
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (int nr) const
2012-06-25 22:20:01 +06:00
{
2015-11-11 22:46:18 +05:00
const Element0d & el = mesh->pointelements[nr];
Ng_Element ret;
ret.type = NG_PNT;
ret.index = el.index;
ret.points.num = 1;
ret.points.ptr = (int*)&el.pnum;
ret.vertices.num = 1;
ret.vertices.ptr = (int*)&el.pnum;
2011-03-09 20:47:04 +05:00
2015-11-11 22:46:18 +05:00
ret.edges.num = 0;
ret.edges.ptr = NULL;
ret.faces.num = 0;
ret.faces.ptr = NULL;
return ret;
}
*/
2014-01-07 16:42:39 +06:00
/*
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
2011-03-09 20:47:04 +05:00
{
const Segment & el = mesh->LineSegment (SegmentIndex(nr));
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&(el[0]);
ret.vertices.num = 2;
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = 1;
ret.edges.ptr = mesh->GetTopology().GetSegmentElementEdgesPtr (nr);
ret.faces.num = 0;
ret.faces.ptr = NULL;
return ret;
}
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (int nr) const
2011-03-09 20:47:04 +05:00
{
const Element2d & el = mesh->SurfaceElement (SurfaceElementIndex (nr));
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&el[0];
ret.vertices.num = el.GetNV();
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = MeshTopology::GetNEdges (el.GetType());
ret.edges.ptr = mesh->GetTopology().GetSurfaceElementEdgesPtr (nr);
ret.faces.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = mesh->GetTopology().GetSurfaceElementFacesPtr (nr);
return ret;
}
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (int nr) const
2011-03-09 20:47:04 +05:00
{
const Element & el = mesh->VolumeElement (ElementIndex (nr));
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&el[0];
ret.vertices.num = el.GetNV();
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = MeshTopology::GetNEdges (el.GetType());
ret.edges.ptr = mesh->GetTopology().GetElementEdgesPtr (nr);
ret.faces.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = mesh->GetTopology().GetElementFacesPtr (nr);
return ret;
}
2014-01-07 16:42:39 +06:00
*/
2011-03-09 20:47:04 +05:00
2015-01-10 19:48:49 +05:00
/*
2013-04-03 02:27:35 +06:00
template <>
DLL_HEADER int Ngx_Mesh :: GetElementIndex<0> (int nr) const
{
return 0;
}
template <>
DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const
{
return (*mesh)[SegmentIndex(nr)].si;
}
template <>
DLL_HEADER int Ngx_Mesh :: GetElementIndex<2> (int nr) const
{
int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex();
return mesh->GetFaceDescriptor(ind).BCProperty();
}
template <>
DLL_HEADER int Ngx_Mesh :: GetElementIndex<3> (int nr) const
{
return (*mesh)[ElementIndex(nr)].GetIndex();
}
2014-01-07 16:42:39 +06:00
*/
2013-04-03 02:27:35 +06:00
/*
2011-03-09 20:47:04 +05:00
DLL_HEADER Ng_Point Ng_GetPoint (int nr)
{
Ng_Point ret;
ret.pt = &mesh->Point(nr + PointIndex::BASE)(0);
return ret;
}
template <>
DLL_HEADER int Ng_GetElementIndex<1> (int nr)
{
return (*mesh)[SegmentIndex(nr)].si;
}
template <>
DLL_HEADER int Ng_GetElementIndex<2> (int nr)
{
int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex();
return mesh->GetFaceDescriptor(ind).BCProperty();
}
template <>
DLL_HEADER int Ng_GetElementIndex<3> (int nr)
{
return (*mesh)[ElementIndex(nr)].GetIndex();
}
2013-04-03 02:27:35 +06:00
template <> int DLL_HEADER Ng_GetNElements<0> ()
{
return 0;
}
template <> int DLL_HEADER Ng_GetNElements<1> ()
{
return mesh->GetNSeg();
}
template <> DLL_HEADER int Ng_GetNElements<2> ()
{
return mesh->GetNSE();
}
template <> DLL_HEADER int Ng_GetNElements<3> ()
{
return mesh->GetNE();
}
template <> DLL_HEADER Ng_Element Ng_GetElement<0> (int nr)
{
cout << "Netgen does not support 0-D elements" << endl;
Ng_Element ret;
return ret;
}
template <> DLL_HEADER Ng_Element Ng_GetElement<1> (int nr)
{
const Segment & el = mesh->LineSegment (SegmentIndex(nr));
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&(el[0]);
ret.vertices.num = 2;
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = 1;
ret.edges.ptr = mesh->GetTopology().GetSegmentElementEdgesPtr (nr);
ret.faces.num = 0;
ret.faces.ptr = NULL;
return ret;
}
template <> DLL_HEADER Ng_Element Ng_GetElement<2> (int nr)
{
const Element2d & el = mesh->SurfaceElement (SurfaceElementIndex (nr));
2011-03-09 20:47:04 +05:00
2013-04-03 02:27:35 +06:00
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&el[0];
ret.vertices.num = el.GetNV();
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = MeshTopology::GetNEdges (el.GetType());
ret.edges.ptr = mesh->GetTopology().GetSurfaceElementEdgesPtr (nr);
ret.faces.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = mesh->GetTopology().GetSurfaceElementFacesPtr (nr);
return ret;
}
template <> DLL_HEADER Ng_Element Ng_GetElement<3> (int nr)
{
const Element & el = mesh->VolumeElement (ElementIndex (nr));
2011-03-09 20:47:04 +05:00
2013-04-03 02:27:35 +06:00
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&el[0];
ret.vertices.num = el.GetNV();
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = MeshTopology::GetNEdges (el.GetType());
ret.edges.ptr = mesh->GetTopology().GetElementEdgesPtr (nr);
ret.faces.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = mesh->GetTopology().GetElementFacesPtr (nr);
return ret;
}
*/
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<3,3> (int elnr,
const double * xi,
double * x,
double * dxdxi) const
{
Point<3> xl(xi[0], xi[1], xi[2]);
Point<3> xg;
Mat<3,3> dx;
mesh->GetCurvedElements().CalcElementTransformation (xl, elnr, xg, dx);
if (x)
for (int i = 0; i < 3; i++) x[i] = xg(i);
if (dxdxi)
for (int i=0; i<3; i++)
{
dxdxi[3*i] = dx(i,0);
dxdxi[3*i+1] = dx(i,1);
dxdxi[3*i+2] = dx(i,2);
}
}
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<2,3> (int elnr,
const double * xi,
double * x,
double * dxdxi) const
{
Point<2> xl(xi[0], xi[1]);
Point<3> xg;
Mat<3,2> dx;
mesh->GetCurvedElements().CalcSurfaceTransformation (xl, elnr, xg, dx);
if (x)
for (int i = 0; i < 3; i++) x[i] = xg(i);
if (dxdxi)
for (int i=0; i<3; i++)
{
dxdxi[2*i] = dx(i,0);
dxdxi[2*i+1] = dx(i,1);
}
}
2016-10-04 22:30:57 +05:00
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<1,3> (int elnr,
const double * xi,
double * x,
double * dxdxi) const
{
Point<3> xg;
Vec<3> dx;
mesh->GetCurvedElements().CalcSegmentTransformation(xi[0],elnr,xg,dx);
if(x)
for(int i=0;i<3;i++) x[i] = xg(i);
if(dxdxi)
for(int i=0;i<3;i++) dxdxi[i] = dx(i);
}
2017-09-22 19:55:10 +05:00
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<0,3> (int elnr,
const double * xi,
double * x,
double * dxdxi) const
{
PointIndex pi = mesh->pointelements[elnr].pnum;
Point<3> xg = mesh->Point(pi);
if (x)
for(int i=0;i<3;i++) x[i] = xg(i);
}
2013-04-03 02:27:35 +06:00
2017-09-22 19:55:10 +05:00
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<2,2> (int elnr,
const double * xi,
double * x,
double * dxdxi) const
{
Point<2> xl(xi[0], xi[1]);
Point<3> xg;
Mat<3,2> dx;
mesh->GetCurvedElements().CalcSurfaceTransformation (xl, elnr, xg, dx);
if (x)
for (int i = 0; i < 2; i++) x[i] = xg(i);
if (dxdxi)
for (int i=0; i<2; i++)
{
dxdxi[2*i] = dx(i,0);
dxdxi[2*i+1] = dx(i,1);
}
}
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<1,2> (int elnr,
const double * xi,
double * x,
double * dxdxi) const
{
Point<3> xg;
Vec<3> dx;
mesh->GetCurvedElements().CalcSegmentTransformation (xi[0], elnr, xg, dx);
if (x)
for (int i = 0; i < 2; i++) x[i] = xg(i);
if (dxdxi)
for (int i=0; i < 2; i++)
dxdxi[i] = dx(i);
}
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<1,1> (int elnr,
const double * xi,
double * x,
double * dxdxi) const
{
2013-05-27 19:01:58 +06:00
Point<3> xg;
Vec<3> dx;
mesh->GetCurvedElements().CalcSegmentTransformation (xi[0], elnr, xg, dx);
if (x) x[0] = xg(0);
if (dxdxi) dxdxi[0] = dx(0);
2013-04-03 02:27:35 +06:00
}
2016-10-16 12:45:16 +05:00
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<0,2> (int elnr,
const double *xi,
double * x,
double * dxdxi) const
{
PointIndex pnum = mesh->pointelements[elnr].pnum;
if (x)
for (int i = 0; i< 2; i++) x[i] = (*mesh)[pnum](i);
}
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER void Ngx_Mesh ::
ElementTransformation<0,1> (int elnr,
const double * xi,
double * x,
double * dxdxi) const
{
2016-08-18 12:50:11 +05:00
PointIndex pnum = mesh->pointelements[elnr].pnum;
2015-11-11 22:46:18 +05:00
if (x) x[0] = (*mesh)[pnum](0);
2016-03-23 23:51:09 +05:00
// if (dxdxi) dxdxi[0] = 0;
// Jacobi-matrix is 1 x 0 !!!
2013-04-03 02:27:35 +06:00
}
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<3,3> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
2011-03-09 20:47:04 +05:00
{
mesh->GetCurvedElements().CalcMultiPointElementTransformation (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
}
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<2,2> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
2011-03-09 20:47:04 +05:00
{
mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<2> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
}
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<2,3> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
2011-03-09 20:47:04 +05:00
{
mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
}
2016-10-16 12:45:16 +05:00
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,3> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
{
mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
}
2017-09-22 19:55:10 +05:00
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,3> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
{
for (int i = 0; i < npts; i++)
ElementTransformation<0,3> (elnr, xi+i*sxi, x+i*sx, dxdxi+i*sdxdxi);
}
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,2> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
2011-03-09 20:47:04 +05:00
{
mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<2> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
}
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,1> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
2011-03-09 20:47:04 +05:00
{
2013-05-27 19:01:58 +06:00
for (int i = 0; i < npts; i++)
ElementTransformation<1,1> (elnr, xi + i*sxi, x+i*sx, dxdxi+i*sdxdxi);
2011-03-09 20:47:04 +05:00
}
2016-10-16 12:45:16 +05:00
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,2> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
{
for (int i = 0; i < npts; i++)
ElementTransformation<0,2> (elnr, xi + i*sxi, x+i*sx, dxdxi+i*sdxdxi);
}
2013-04-03 02:27:35 +06:00
template <> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,1> (int elnr, int npts,
const double * xi, size_t sxi,
double * x, size_t sx,
double * dxdxi, size_t sdxdxi) const
2012-07-09 15:50:48 +06:00
{
2015-11-11 22:46:18 +05:00
for (int i = 0; i < npts; i++)
ElementTransformation<0,1> (elnr, xi + i*sxi, x+i*sx, dxdxi+i*sdxdxi);
2012-07-09 15:50:48 +06:00
}
2011-03-09 20:47:04 +05:00
void Ngx_Mesh :: GetParentNodes (int ni, int * parents) const
{
ni++;
if (ni <= mesh->mlbetweennodes.Size())
{
parents[0] = mesh->mlbetweennodes.Get(ni).I1()-1;
parents[1] = mesh->mlbetweennodes.Get(ni).I2()-1;
}
else
parents[0] = parents[1] = -1;
}
int Ngx_Mesh :: GetParentElement (int ei) const
{
ei++;
if (mesh->GetDimension() == 3)
{
if (ei <= mesh->mlparentelement.Size())
return mesh->mlparentelement.Get(ei)-1;
}
else
{
if (ei <= mesh->mlparentsurfaceelement.Size())
return mesh->mlparentsurfaceelement.Get(ei)-1;
}
return -1;
}
int Ngx_Mesh :: GetParentSElement (int ei) const
{
ei++;
if (mesh->GetDimension() == 3)
{
if (ei <= mesh->mlparentsurfaceelement.Size())
return mesh->mlparentsurfaceelement.Get(ei)-1;
}
else
{
return -1;
}
return -1;
}
2011-03-09 20:47:04 +05:00
2013-04-03 02:27:35 +06:00
2016-04-10 09:36:05 +05:00
2017-11-27 14:19:43 +05:00
#ifdef __SSE__
#include <immintrin.h>
2016-04-10 09:36:05 +05:00
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,1> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd * xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2016-04-10 09:36:05 +05:00
{
cout << "multi-eltrafo simd called, 1,1,simd" << endl;
}
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<2,2> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd * xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2016-04-10 09:36:05 +05:00
{
2016-07-11 21:27:44 +05:00
mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<2>
(elnr, npts,
reinterpret_cast<const SIMD<double>*> (xi), sxi,
reinterpret_cast<SIMD<double>*> (x), sx,
reinterpret_cast<SIMD<double>*> (dxdxi), sdxdxi);
2016-07-11 21:27:44 +05:00
/*
2016-04-10 09:36:05 +05:00
for (int i = 0; i < npts; i++)
{
double hxi[4][2];
double hx[4][2];
double hdxdxi[4][4];
for (int j = 0; j < 4; j++)
for (int k = 0; k < 2; k++)
hxi[j][k] = ((double*)&(xi[k]))[j];
MultiElementTransformation<2,2> (elnr, 4, &hxi[0][0], 2, &hx[0][0], 2, &hdxdxi[0][0], 4);
for (int j = 0; j < 4; j++)
for (int k = 0; k < 2; k++)
((double*)&(x[k]))[j] = hx[j][k];
for (int j = 0; j < 4; j++)
for (int k = 0; k < 4; k++)
((double*)&(dxdxi[k]))[j] = hdxdxi[j][k];
xi += sxi;
x += sx;
dxdxi += sdxdxi;
}
2016-07-11 21:27:44 +05:00
*/
2016-04-10 09:36:05 +05:00
}
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<3,3> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd * xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2016-04-10 09:36:05 +05:00
{
mesh->GetCurvedElements().CalcMultiPointElementTransformation
(elnr, npts,
reinterpret_cast<const SIMD<double>*> (xi), sxi,
reinterpret_cast<SIMD<double>*> (x), sx,
reinterpret_cast<SIMD<double>*> (dxdxi), sdxdxi);
/*
2016-04-10 09:36:05 +05:00
for (int i = 0; i < npts; i++)
{
double hxi[4][3];
double hx[4][3];
double hdxdxi[4][9];
for (int j = 0; j < 4; j++)
for (int k = 0; k < 3; k++)
hxi[j][k] = ((double*)&(xi[k]))[j];
MultiElementTransformation<3,3> (elnr, 4, &hxi[0][0], 3, &hx[0][0], 3, &hdxdxi[0][0], 9);
for (int j = 0; j < 4; j++)
for (int k = 0; k < 3; k++)
((double*)&(x[k]))[j] = hx[j][k];
for (int j = 0; j < 4; j++)
for (int k = 0; k < 9; k++)
((double*)&(dxdxi[k]))[j] = hdxdxi[j][k];
xi += sxi;
x += sx;
dxdxi += sdxdxi;
}
*/
2016-04-10 09:36:05 +05:00
}
2016-10-16 12:45:16 +05:00
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,2> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd *xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2016-10-16 12:45:16 +05:00
{
cout << "MultiElementtransformation<0,2> simd not implemented" << endl;
}
2016-04-10 09:36:05 +05:00
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,1> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd * xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2016-04-10 09:36:05 +05:00
{
cout << "multi-eltrafo simd called, 0,1,simd" << endl;
}
2016-10-16 12:45:16 +05:00
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,3> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd * xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2016-10-16 12:45:16 +05:00
{
2017-11-24 01:26:23 +05:00
mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<3>
(elnr, npts,
reinterpret_cast<const SIMD<double>*> (xi), sxi,
reinterpret_cast<SIMD<double>*> (x), sx,
reinterpret_cast<SIMD<double>*> (dxdxi), sdxdxi);
/*
2016-10-16 12:45:16 +05:00
double hxi[4][1];
double hx[4][3];
double hdxdxi[4][3];
for (int j = 0; j<4;j++)
hxi[j][0] = ((double*)&(xi[0]))[j];
MultiElementTransformation<1,3> (elnr, 4, &hxi[0][0], 1, &hx[0][0], 3, &hdxdxi[0][0],3);
2016-10-16 12:45:16 +05:00
for(int j=0; j<4; j++)
for(int k=0; k<3; k++)
((double*)&(x[k]))[j] = hx[j][k];
for(int j=0; j< 4; j++)
for (int k = 0; k<3; k++)
((double*) & (dxdxi[k]))[j] = hdxdxi[j][k];
xi += sxi;
x += sx;
dxdxi += sdxdxi;
2017-11-24 01:26:23 +05:00
*/
2016-10-16 12:45:16 +05:00
}
2016-04-10 09:36:05 +05:00
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,2> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd * xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2016-04-10 09:36:05 +05:00
{
2017-11-24 01:26:23 +05:00
mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<2>
(elnr, npts,
reinterpret_cast<const SIMD<double>*> (xi), sxi,
reinterpret_cast<SIMD<double>*> (x), sx,
reinterpret_cast<SIMD<double>*> (dxdxi), sdxdxi);
/*
2016-05-17 18:35:24 +05:00
for (int i = 0; i < npts; i++)
{
double hxi[4][1];
double hx[4][2];
2017-04-12 19:44:02 +05:00
double hdxdxi[4][2];
2016-05-17 18:35:24 +05:00
for (int j = 0; j < 4; j++)
for (int k = 0; k < 1; k++)
hxi[j][k] = ((double*)&(xi[k]))[j];
2017-04-12 19:44:02 +05:00
MultiElementTransformation<1,2> (elnr, 4, &hxi[0][0], 1, &hx[0][0], 2, &hdxdxi[0][0], 2);
2016-05-17 18:35:24 +05:00
for (int j = 0; j < 4; j++)
for (int k = 0; k < 2; k++)
((double*)&(x[k]))[j] = hx[j][k];
for (int j = 0; j < 4; j++)
2017-04-12 19:44:02 +05:00
for (int k = 0; k < 2; k++)
2016-05-17 18:35:24 +05:00
((double*)&(dxdxi[k]))[j] = hdxdxi[j][k];
xi += sxi;
x += sx;
dxdxi += sdxdxi;
}
2017-11-24 01:26:23 +05:00
*/
2016-04-10 09:36:05 +05:00
}
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<2,3> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd * xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2016-04-10 09:36:05 +05:00
{
2016-07-12 15:45:21 +05:00
mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3>
(elnr, npts,
reinterpret_cast<const SIMD<double>*> (xi), sxi,
reinterpret_cast<SIMD<double>*> (x), sx,
reinterpret_cast<SIMD<double>*> (dxdxi), sdxdxi);
2016-07-12 15:45:21 +05:00
/*
2016-05-27 21:06:20 +05:00
for (int i = 0; i < npts; i++)
{
double hxi[4][2];
double hx[4][3];
double hdxdxi[4][6];
for (int j = 0; j < 4; j++)
for (int k = 0; k < 2; k++)
hxi[j][k] = ((double*)&(xi[k]))[j];
MultiElementTransformation<2,3> (elnr, 4, &hxi[0][0], 2, &hx[0][0], 3, &hdxdxi[0][0], 6);
for (int j = 0; j < 4; j++)
for (int k = 0; k < 3; k++)
((double*)&(x[k]))[j] = hx[j][k];
for (int j = 0; j < 4; j++)
for (int k = 0; k < 6; k++)
((double*)&(dxdxi[k]))[j] = hdxdxi[j][k];
xi += sxi;
x += sx;
dxdxi += sdxdxi;
}
2016-07-12 15:45:21 +05:00
*/
2016-04-10 09:36:05 +05:00
}
2017-09-22 19:55:10 +05:00
template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,3> (int elnr, int npts,
2017-11-24 01:26:23 +05:00
const tAVXd * xi, size_t sxi,
tAVXd * x, size_t sx,
tAVXd * dxdxi, size_t sdxdxi) const
2017-09-22 19:55:10 +05:00
{
for (int i = 0; i < npts; i++)
{
double hxi[4][1];
double hx[4][3];
for (int j = 0; j < 4; j++)
for (int k = 0; k < 1; k++)
hxi[j][k] = ((double*)&(xi[k]))[j];
MultiElementTransformation<0,3> (elnr, 4, &hxi[0][0], 2, &hx[0][0], 3, (double*)nullptr, 0);
for (int j = 0; j < 4; j++)
for (int k = 0; k < 3; k++)
((double*)&(x[k]))[j] = hx[j][k];
xi += sxi;
x += sx;
dxdxi += sdxdxi;
}
}
#endif
2016-04-10 09:36:05 +05:00
2015-08-31 20:41:26 +05:00
template <>
DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <1>
(double * hp, double * lami,
bool build_searchtree,
int * const indices, int numind) const
{
if (mesh->GetDimension() != 1)
throw NgException("FindElementOfPoint<1> called for multidim mesh");
Point<3> p(hp[0], 0,0);
for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
{
auto & seg = (*mesh)[si];
Point<3> p1 = (*mesh)[seg[0]];
Point<3> p2 = (*mesh)[seg[1]];
double lam = (p(0)-p1(0)) / (p2(0)-p1(0));
if (lam >= -1e-10 && lam <= 1+1e-10)
{
lami[0] = 1-lam;
return si;
}
}
return -1;
}
2013-04-03 02:27:35 +06:00
template <>
2014-01-07 16:42:39 +06:00
DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <2>
2013-04-03 02:27:35 +06:00
(double * p, double * lami,
bool build_searchtree,
2014-01-07 16:42:39 +06:00
int * const indices, int numind) const
2013-04-03 02:27:35 +06:00
{
Array<int> dummy(numind);
for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
double lam3[3];
2016-05-09 08:45:09 +05:00
int ind;
if (mesh->GetDimension() == 2)
{
Point<3> p2d(p[0], p[1], 0);
ind = mesh->GetElementOfPoint(p2d, lam3, &dummy, build_searchtree);
}
else
{
Point3d p3d(p[0], p[1], p[2]);
ind = mesh->GetSurfaceElementOfPoint(p3d, lam3, &dummy, build_searchtree);
}
2013-04-03 02:27:35 +06:00
if (ind > 0)
{
if(mesh->SurfaceElement(ind).GetType()==QUAD)
{
lami[0] = lam3[0];
lami[1] = lam3[1];
}
else
{
lami[0] = 1-lam3[0]-lam3[1];
lami[1] = lam3[0];
}
}
return ind-1;
}
template <>
DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <3>
(double * p, double * lami,
bool build_searchtree,
2014-01-07 16:42:39 +06:00
int * const indices, int numind) const
2013-04-03 02:27:35 +06:00
{
Array<int> dummy(numind);
for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
Point<3> p3d(p[0], p[1], p[2]);
int ind =
mesh->GetElementOfPoint(p3d, lami, &dummy, build_searchtree);
return ind-1;
}
2018-05-02 00:20:54 +05:00
void Ngx_Mesh :: Curve (int order)
{
NgLock meshlock (mesh->MajorMutex(), true);
mesh->BuildCurvedElements(order);
}
void Ngx_Mesh :: Refine (NG_REFINEMENT_TYPE reftype,
2018-01-04 14:43:22 +05:00
void (*task_manager)(function<void(int,int)>),
Tracer tracer)
{
NgLock meshlock (mesh->MajorMutex(), 1);
BisectionOptions biopt;
biopt.usemarkedelements = 1;
biopt.refine_p = 0;
biopt.refine_hp = 0;
if (reftype == NG_REFINE_P)
biopt.refine_p = 1;
if (reftype == NG_REFINE_HP)
biopt.refine_hp = 1;
biopt.task_manager = task_manager;
2018-01-04 14:43:22 +05:00
biopt.tracer = tracer;
const Refinement & ref = mesh->GetGeometry()->GetRefinement();
ref.Bisect (*mesh, biopt);
2018-01-04 14:43:22 +05:00
(*tracer)("call updatetop", false);
2018-01-04 22:45:07 +05:00
mesh -> UpdateTopology(task_manager, tracer);
2018-01-04 14:43:22 +05:00
(*tracer)("call updatetop", true);
mesh -> GetCurvedElements().SetIsHighOrder (false);
}
2013-04-03 02:27:35 +06:00
2014-12-03 02:53:29 +05:00
#ifdef PARALLEL
2014-12-03 02:50:38 +05:00
std::tuple<int,int*> Ngx_Mesh :: GetDistantProcs (int nodetype, int locnum) const
{
switch (nodetype)
{
case 0:
{
FlatArray<int> dn = mesh->GetParallelTopology().GetDistantPNums(locnum);
return std::tuple<int,int*>(dn.Size(), &dn[0]);
}
case 1:
{
FlatArray<int> dn = mesh->GetParallelTopology().GetDistantEdgeNums(locnum);
return std::tuple<int,int*>(dn.Size(), &dn[0]);
}
case 2:
{
FlatArray<int> dn = mesh->GetParallelTopology().GetDistantFaceNums(locnum);
return std::tuple<int,int*>(dn.Size(), &dn[0]);
}
default:
return std::tuple<int,int*>(0,nullptr);
}
}
2013-04-03 02:27:35 +06:00
2014-12-03 02:53:29 +05:00
#endif
2011-03-09 20:47:04 +05:00
}
int link_it_nginterface_v2;