Commit 60f42ec0 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/surdDist-general' into 'master'

Extend SurfDistance module of FieldConvert to Quad and Hex elements

This MR extends the SurfDistance module of FieldConvert to quadrilateral and hexahedron elements.

I tested it with quads and with a simple hex mesh. It would be nice if someone could test it with a mesh with a prism boundary layer.

See merge request !524
parents f8efa57e 4195eaf7
......@@ -19,6 +19,7 @@ v4.3.0
- Fix restarting from checkpoint file (!517)
**FieldConvert:**
- Extended surface distance module to support hexahedra and quads (!524)
- Small fixes in interpolation routine (!515)
- Add support for surface extraction in 3DH1D case (!521)
- Add support for isocontour extraction for 3DH1D (!525)
......
......@@ -83,6 +83,7 @@ ADD_NEKTAR_TEST(bfs_vort)
ADD_NEKTAR_TEST(bfs_vort_rng)
# ADD_NEKTAR_TEST(chan3D_pts)
ADD_NEKTAR_TEST(chan3DH1D_meanmode)
ADD_NEKTAR_TEST(cube_prismhex)
# windows produces slightly differently formatted files which results in
# different hashes
......
......@@ -67,6 +67,7 @@ void ProcessSurfDistance::Process(po::variables_map &vm)
{
int i, j, k, cnt;
int surf = m_config["bnd"].as<int>();
int expdim = m_f->m_graph->GetMeshDimension();
ASSERTL0(surf >= 0, "Invalid surface "+boost::lexical_cast<string>(surf));
......@@ -95,6 +96,9 @@ void ProcessSurfDistance::Process(po::variables_map &vm)
m_f->m_fielddef[0]->m_fields[0] = "dist";
}
ASSERTL0(!(m_f->m_fielddef[0]->m_numHomogeneousDir),
"Homogeneous expansions not supported");
for (i = cnt = 0; i < BndExp.num_elements(); ++i)
{
if (i != surf)
......@@ -107,17 +111,70 @@ void ProcessSurfDistance::Process(po::variables_map &vm)
{
int elmtNum = BoundarytoElmtID [cnt];
int facetNum = BoundarytoTraceID[cnt];
int oppositeNum;
// Get boundary and element expansions.
LocalRegions::ExpansionSharedPtr bndElmt = BndExp[i]->GetExp(j);
LocalRegions::ExpansionSharedPtr elmt =
m_f->m_exp[0]->GetExp(elmtNum);
ASSERTL0(elmt->DetShapeType() == LibUtilities::ePrism,
"Only prisms supported for now!");
// Determine which face is opposite to the surface
switch(elmt->DetShapeType())
{
case LibUtilities::eQuadrilateral:
{
oppositeNum = (facetNum+2)%4;
}
break;
ASSERTL0(facetNum == 1 || facetNum == 3,
case LibUtilities::ePrism:
{
switch(facetNum)
{
case 1:
oppositeNum = 3;
break;
case 3:
oppositeNum = 1;
break;
default:
ASSERTL0(false,
"Surface must be on a triangular face of the prism.");
}
}
break;
case LibUtilities::eHexahedron:
{
switch(facetNum)
{
case 0:
oppositeNum = 5;
break;
case 1:
oppositeNum = 3;
break;
case 2:
oppositeNum = 4;
break;
case 3:
oppositeNum = 1;
break;
case 4:
oppositeNum = 2;
break;
case 5:
oppositeNum = 0;
break;
default:
ASSERTL0(false, "Face out of bound.");
}
}
break;
default:
ASSERTL0(false, "Element not supported");
}
int nq = elmt ->GetTotPoints();
int nqBnd = bndElmt->GetTotPoints();
......@@ -128,21 +185,48 @@ void ProcessSurfDistance::Process(po::variables_map &vm)
x[2] = Array<OneD, NekDouble>(nq);
elmt->GetCoords(x[0], x[1], x[2]);
Array<OneD, NekDouble> face1(nqBnd), face3(nqBnd);
Array<OneD, NekDouble> face1(nqBnd), face2(nqBnd);
Array<OneD, NekDouble> dist =
BndExp[i]->UpdatePhys() + BndExp[i]->GetPhys_Offset(j);
// Zero existing value.
Vmath::Zero(nqBnd, dist, 1);
// Calculate distance between two faces of prism.
for (k = 0; k < 3; ++k)
// Calculate distance between two faces of the element
for (k = 0; k < expdim; ++k)
{
switch(expdim)
{
case 2:
{
elmt->GetEdgePhysVals(facetNum, bndElmt, x[k], face1);
elmt->GetEdgePhysVals(oppositeNum, bndElmt, x[k], face2);
// Consider edge orientation
if (elmt->GetEorient(facetNum) ==
elmt->GetEorient(oppositeNum))
{
elmt->GetFacePhysVals(1, bndElmt, x[k], face1);
elmt->GetFacePhysVals(3, bndElmt, x[k], face3);
Vmath::Vsub (nqBnd, face1, 1, face3, 1, face1, 1);
Vmath::Reverse(nqBnd, face2, 1, face2, 1);
}
}
break;
case 3:
{
// Use orientation from the surface for both faces
StdRegions::Orientation orientation =
elmt->GetForient(facetNum);
elmt->GetFacePhysVals(facetNum, bndElmt,
x[k], face1, orientation);
elmt->GetFacePhysVals(oppositeNum, bndElmt,
x[k], face2, orientation);
}
break;
default:
ASSERTL0(false, "Expansion not supported");
}
Vmath::Vsub (nqBnd, face1, 1, face2, 1, face1, 1);
Vmath::Vvtvp(nqBnd, face1, 1, face1, 1, dist, 1, dist, 1);
}
Vmath::Vsqrt(nqBnd, dist, 1, dist, 1);
}
BndExp[i]->FwdTrans(BndExp[i]->GetPhys(), BndExp[i]->UpdateCoeffs());
......
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Surface distance calculation on 3D hex/prism mesh</description>
<executable>FieldConvert</executable>
<parameters>-e -m surfdistance:bnd=0 cube_prismhex.xml out.fld</parameters>
<files>
<file description="Session File">cube_prismhex.xml</file>
</files>
<metrics>
<metric type="l2" id="1">
<value variable="dist" tolerance="1e-12">0.5</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR>
<GEOMETRY DIM="3" SPACE="3">
<VERTEX>
<V ID="0">0.00000000e+00 5.00000000e-01 0.00000000e+00</V>
<V ID="1">5.00000000e-01 5.00000000e-01 0.00000000e+00</V>
<V ID="2">5.00000000e-01 1.00000000e+00 0.00000000e+00</V>
<V ID="3">0.00000000e+00 1.00000000e+00 0.00000000e+00</V>
<V ID="4">0.00000000e+00 5.00000000e-01 5.00000000e-01</V>
<V ID="5">5.00000000e-01 5.00000000e-01 5.00000000e-01</V>
<V ID="6">5.00000000e-01 1.00000000e+00 5.00000000e-01</V>
<V ID="7">0.00000000e+00 1.00000000e+00 5.00000000e-01</V>
<V ID="8">0.00000000e+00 5.00000000e-01 1.00000000e+00</V>
<V ID="9">5.00000000e-01 5.00000000e-01 1.00000000e+00</V>
<V ID="10">5.00000000e-01 1.00000000e+00 1.00000000e+00</V>
<V ID="11">0.00000000e+00 1.00000000e+00 1.00000000e+00</V>
<V ID="12">1.00000000e+00 5.00000000e-01 0.00000000e+00</V>
<V ID="13">1.00000000e+00 1.00000000e+00 0.00000000e+00</V>
<V ID="14">1.00000000e+00 5.00000000e-01 5.00000000e-01</V>
<V ID="15">1.00000000e+00 1.00000000e+00 5.00000000e-01</V>
<V ID="16">1.00000000e+00 5.00000000e-01 1.00000000e+00</V>
<V ID="17">1.00000000e+00 1.00000000e+00 1.00000000e+00</V>
<V ID="18">0.00000000e+00 0.00000000e+00 5.00000000e-01</V>
<V ID="19">5.00000000e-01 0.00000000e+00 5.00000000e-01</V>
<V ID="20">5.00000000e-01 0.00000000e+00 0.00000000e+00</V>
<V ID="21">0.00000000e+00 0.00000000e+00 0.00000000e+00</V>
<V ID="22">0.00000000e+00 0.00000000e+00 1.00000000e+00</V>
<V ID="23">5.00000000e-01 0.00000000e+00 1.00000000e+00</V>
<V ID="24">1.00000000e+00 0.00000000e+00 5.00000000e-01</V>
<V ID="25">1.00000000e+00 0.00000000e+00 0.00000000e+00</V>
<V ID="26">1.00000000e+00 0.00000000e+00 1.00000000e+00</V>
</VERTEX>
<EDGE>
<E ID="0"> 0 1 </E>
<E ID="1"> 1 2 </E>
<E ID="2"> 2 3 </E>
<E ID="3"> 3 0 </E>
<E ID="4"> 0 4 </E>
<E ID="5"> 1 5 </E>
<E ID="6"> 2 6 </E>
<E ID="7"> 3 7 </E>
<E ID="8"> 4 5 </E>
<E ID="9"> 5 6 </E>
<E ID="10"> 6 7 </E>
<E ID="11"> 7 4 </E>
<E ID="12"> 4 8 </E>
<E ID="13"> 5 9 </E>
<E ID="14"> 6 10 </E>
<E ID="15"> 7 11 </E>
<E ID="16"> 8 9 </E>
<E ID="17"> 9 10 </E>
<E ID="18"> 10 11 </E>
<E ID="19"> 11 8 </E>
<E ID="20"> 1 12 </E>
<E ID="21"> 12 13 </E>
<E ID="22"> 13 2 </E>
<E ID="23"> 12 14 </E>
<E ID="24"> 13 15 </E>
<E ID="25"> 5 14 </E>
<E ID="26"> 14 15 </E>
<E ID="27"> 15 6 </E>
<E ID="28"> 14 16 </E>
<E ID="29"> 15 17 </E>
<E ID="30"> 9 16 </E>
<E ID="31"> 16 17 </E>
<E ID="32"> 17 10 </E>
<E ID="33"> 18 19 </E>
<E ID="34"> 19 20 </E>
<E ID="35"> 21 20 </E>
<E ID="36"> 18 21 </E>
<E ID="37"> 18 5 </E>
<E ID="38"> 19 5 </E>
<E ID="39"> 20 1 </E>
<E ID="40"> 21 1 </E>
<E ID="41"> 22 23 </E>
<E ID="42"> 23 19 </E>
<E ID="43"> 22 18 </E>
<E ID="44"> 22 9 </E>
<E ID="45"> 23 9 </E>
<E ID="46"> 18 4 </E>
<E ID="47"> 21 0 </E>
<E ID="48"> 22 8 </E>
<E ID="49"> 24 14 </E>
<E ID="50"> 25 12 </E>
<E ID="51"> 24 25 </E>
<E ID="52"> 19 24 </E>
<E ID="53"> 19 14 </E>
<E ID="54"> 20 12 </E>
<E ID="55"> 20 25 </E>
<E ID="56"> 26 16 </E>
<E ID="57"> 26 24 </E>
<E ID="58"> 23 26 </E>
<E ID="59"> 23 16 </E>
</EDGE>
<FACE>
<Q ID="0"> 0 1 2 3</Q>
<Q ID="1"> 0 5 8 4</Q>
<Q ID="2"> 1 6 9 5</Q>
<Q ID="3"> 2 6 10 7</Q>
<Q ID="4"> 3 7 11 4</Q>
<Q ID="5"> 8 9 10 11</Q>
<Q ID="6"> 8 13 16 12</Q>
<Q ID="7"> 9 14 17 13</Q>
<Q ID="8"> 10 14 18 15</Q>
<Q ID="9"> 11 15 19 12</Q>
<Q ID="10"> 16 17 18 19</Q>
<Q ID="11"> 20 21 22 1</Q>
<Q ID="12"> 20 23 25 5</Q>
<Q ID="13"> 21 24 26 23</Q>
<Q ID="14"> 22 24 27 6</Q>
<Q ID="15"> 25 26 27 9</Q>
<Q ID="16"> 25 28 30 13</Q>
<Q ID="17"> 26 29 31 28</Q>
<Q ID="18"> 27 29 32 14</Q>
<Q ID="19"> 30 31 32 17</Q>
<Q ID="20"> 33 34 35 36</Q>
<T ID="21"> 33 38 37</T>
<Q ID="22"> 34 39 5 38</Q>
<T ID="23"> 35 39 40</T>
<Q ID="24"> 36 40 5 37</Q>
<Q ID="25"> 41 42 33 43</Q>
<T ID="26"> 41 45 44</T>
<Q ID="27"> 42 38 13 45</Q>
<Q ID="28"> 43 37 13 44</Q>
<Q ID="29"> 46 36 47 4</Q>
<T ID="30"> 46 37 8</T>
<T ID="31"> 47 40 0</T>
<Q ID="32"> 48 43 46 12</Q>
<T ID="33"> 48 44 16</T>
<Q ID="34"> 49 23 50 51</Q>
<T ID="35"> 49 53 52</T>
<Q ID="36"> 23 54 34 53</Q>
<T ID="37"> 50 54 55</T>
<Q ID="38"> 51 55 34 52</Q>
<Q ID="39"> 56 28 49 57</Q>
<T ID="40"> 56 59 58</T>
<Q ID="41"> 28 53 42 59</Q>
<Q ID="42"> 57 52 42 58</Q>
<T ID="43"> 53 25 38</T>
<T ID="44"> 54 20 39</T>
<T ID="45"> 59 30 45</T>
</FACE>
<ELEMENT>
<H ID="0"> 0 1 2 3 4 5 </H>
<H ID="1"> 5 6 7 8 9 10 </H>
<H ID="2"> 11 12 13 14 2 15 </H>
<H ID="3"> 15 16 17 18 7 19 </H>
<R ID="4"> 20 21 22 23 24 </R>
<R ID="5"> 25 26 27 21 28 </R>
<R ID="6"> 29 30 24 31 1 </R>
<R ID="7"> 32 33 28 30 6 </R>
<R ID="8"> 34 35 36 37 38 </R>
<R ID="9"> 39 40 41 35 42 </R>
<R ID="10"> 36 43 12 44 22 </R>
<R ID="11"> 41 45 16 43 27 </R>
</ELEMENT>
<COMPOSITE>
<C ID="0"> H[0-3] </C>
<C ID="1"> F[26,33,40,45,10,19] </C>
<C ID="2"> R[4-11] </C>
</COMPOSITE>
<DOMAIN> C[0,2] </DOMAIN>
</GEOMETRY>
<EXPANSIONS>
<E COMPOSITE="C[0]" NUMMODES="4" TYPE="MODIFIED" FIELDS="u" />
<E COMPOSITE="C[2]" NUMMODES="4" TYPE="MODIFIED" FIELDS="u" />
</EXPANSIONS>
<CONDITIONS>
<VARIABLES>
<V ID="0"> u </V>
</VARIABLES>
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" VALUE="0" />
</REGION>
</BOUNDARYCONDITIONS>
</CONDITIONS>
</NEKTAR>
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment