mesh.SecondOrder : fix Segment mapping

This commit is contained in:
Joachim Schoeberl 2021-10-15 18:52:11 +02:00
parent 04a31eeed3
commit 95301e11ba
4 changed files with 36 additions and 1 deletions

View File

@ -43,7 +43,7 @@ namespace netgen
Point (const Point<D,T2> & p2) Point (const Point<D,T2> & p2)
{ for (int i = 0; i < D; i++) x[i] = p2(i); } { for (int i = 0; i < D; i++) x[i] = p2(i); }
explicit Point (const Vec<D> & v) explicit Point (const Vec<D,T> & v)
{ for (int i = 0; i < D; i++) x[i] = v(i); } { for (int i = 0; i < D; i++) x[i] = v(i); }

View File

@ -78,6 +78,15 @@ namespace netgen
return res; return res;
} }
template <int D, int SW>
inline Vec<D, SIMD<double,SW>> operator* (SIMD<double,SW> s, Vec<D,double> b)
{
Vec<D,SIMD<double,SW>> res;
for (int i = 0; i < D; i++)
res(i) = s * b(i);
return res;
}
template <int D> template <int D>
inline double operator* (Vec<D> a, Vec<D> b) inline double operator* (Vec<D> a, Vec<D> b)

View File

@ -1429,6 +1429,31 @@ namespace netgen
info.order = order; info.order = order;
info.ndof = info.nv = 2; info.ndof = info.nv = 2;
if (order == 1)
{
auto & seg = mesh.LineSegment(elnr);
if (seg.GetNP() == 2)
{
if (x)
*x = Point<3,T>((1-xi) * Vec<3>(mesh[seg.PNums()[0]]) + xi * Vec<3>(mesh[seg.PNums()[1]]));
if (dxdxi)
*dxdxi = Vec<3>(mesh[seg.PNums()[1]])-Vec<3>(mesh[seg.PNums()[0]]);
}
else
{
if (x)
{
*x = Point<3,T>(2*(xi-1)*(xi-0.5) * Vec<3>(mesh[seg.PNums()[1]])
+ 4*xi*(1-xi) * Vec<3>(mesh[seg.PNums()[2]])
+ 2*xi*(xi-0.5) * Vec<3>(mesh[seg.PNums()[0]]));
}
if (dxdxi)
*dxdxi = T(1.0) * (Vec<3>(mesh[seg.PNums()[1]])-Vec<3>(mesh[seg.PNums()[0]]))
+ (4-8*xi)*Vec<3>(mesh[seg.PNums()[2]]);
}
return;
}
if (info.order > 1) if (info.order > 1)
{ {
const MeshTopology & top = mesh.GetTopology(); const MeshTopology & top = mesh.GetTopology();

View File

@ -110,6 +110,7 @@ namespace netgen
EDGEPOINT); EDGEPOINT);
between.Set (i2, el[2]); between.Set (i2, el[2]);
} }
el.SetCurved(true);
} }
// refine surface elements // refine surface elements