Commit a06e0c79 authored by Michael Turner's avatar Michael Turner

remove volume curvature work

parent b52ac3b0
......@@ -33,12 +33,8 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/Foundations/Interp.h>
#include <StdRegions/StdNodalTetExp.h>
#include <SpatialDomains/Geometry3D.h>
#include <SpatialDomains/Geometry2D.h>
#include <SpatialDomains/GeomFactors.h>
#include <SpatialDomains/PointGeom.h>
......@@ -255,76 +251,6 @@ namespace SpatialDomains
int i,j,k;
if (m_curve)
{
int pdim = LibUtilities::PointsManager()[
LibUtilities::PointsKey(3, m_curve->m_ptype)]
->GetPointsDim();
// Deal with 2D points type separately
// (e.g. electrostatic or Fekete points) to 1D tensor
// product.
if (pdim == 3)
{
int N = m_curve->m_points.size();
// Need the cubic formula to solve for nEdgePts
NekDouble tmp1 = pow(sqrt(3.0*(243.0*N*N - 1.0)) + 27.0*N, 1.0/3.0);
int nEdgePts = round(tmp1 / pow(3.0, 2.0/3.0) +
1.0 / tmp1 / pow(3.0, 1.0/3.0) - 2.0) + 1;
ASSERTL0(nEdgePts*(nEdgePts+1)*(nEdgePts+2) / 6 == N,
"NUMPOINTS should be a tetrahedral number for"
" tetrahedron "
+ boost::lexical_cast<string>(m_globalID));
const LibUtilities::PointsKey P0(
nEdgePts, LibUtilities::eGaussLobattoLegendre);
const LibUtilities::PointsKey P1(
nEdgePts, LibUtilities::eGaussRadauMAlpha1Beta0);
const LibUtilities::PointsKey P2(
nEdgePts, LibUtilities::eGaussRadauMAlpha2Beta0);
const LibUtilities::BasisKey T0(
LibUtilities::eOrtho_A, nEdgePts, P0);
const LibUtilities::BasisKey T1(
LibUtilities::eOrtho_B, nEdgePts, P1);
const LibUtilities::BasisKey T2(
LibUtilities::eOrtho_C, nEdgePts, P2);
Array<OneD, NekDouble> phys(max(N, m_xmap->GetTotPoints()));
Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
for(i = 0; i < m_coordim; ++i)
{
// Create a StdNodalTetExp.
StdRegions::StdNodalTetExpSharedPtr t =
MemoryManager<StdRegions::StdNodalTetExp>
::AllocateSharedPtr(T0, T1, T2, m_curve->m_ptype);
for (j = 0; j < N; ++j)
{
phys[j] = (m_curve->m_points[j]->GetPtr())[i];
}
t->BwdTrans(phys, tmp);
// Interpolate points to standard region.
LibUtilities::Interp3D(
P0, P1, P2, tmp,
m_xmap->GetBasis(0)->GetPointsKey(),
m_xmap->GetBasis(1)->GetPointsKey(),
m_xmap->GetBasis(2)->GetPointsKey(),
phys);
// Forwards transform to get coefficient space.
m_xmap->FwdTrans(phys, m_coeffs[i]);
}
}
else
{
ASSERTL0(false, "Only 3D points distributions supported.");
}
}
for(i = 0; i < m_forient.size(); i++)
{
m_faces[i]->FillGeom();
......
......@@ -95,7 +95,6 @@ namespace Nektar
bool m_owndata;
int m_eid;
bool m_ownverts;
CurveSharedPtr m_curve;
//---------------------------------------
// 3D Geometry Methods
......
This diff is collapsed.
......@@ -390,7 +390,6 @@ namespace Nektar
SPATIAL_DOMAINS_EXPORT CurveMap& GetCurvedEdges() { return m_curvedEdges; }
SPATIAL_DOMAINS_EXPORT CurveMap& GetCurvedFaces() { return m_curvedFaces; }
SPATIAL_DOMAINS_EXPORT CurveMap& GetCurvedVolumes() { return m_curvedVolumes; }
// void AddExpansion(ExpansionShPtr expansion) { m_expansions[expansion->m_geomShPtr->GetGlobalID()] = expansion; }
SPATIAL_DOMAINS_EXPORT const PointGeomMap& GetAllPointGeoms() const { return m_vertSet; }
......@@ -413,7 +412,6 @@ namespace Nektar
CurveMap m_curvedEdges;
CurveMap m_curvedFaces;
CurveMap m_curvedVolumes;
SegGeomMap m_segGeoms;
......
......@@ -89,7 +89,7 @@ namespace Nektar
mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
ReadCurves(doc);
ReadEdges(doc);
ReadElements(doc);
......@@ -111,9 +111,9 @@ namespace Nektar
ASSERTL0(field, "Unable to find EDGE tag in file.");
string IsCompressed;
field->QueryStringAttribute("COMPRESSED",&IsCompressed);
if(IsCompressed.size())
field->QueryStringAttribute("COMPRESSED",&IsCompressed);
if(IsCompressed.size())
{
ASSERTL0(boost::iequals(IsCompressed,
LibUtilities::CompressData::GetCompressString()),
......@@ -165,7 +165,7 @@ namespace Nektar
/// ? being the element type.
/// Read the ID field first.
TiXmlElement *edge = field->FirstChildElement("E");
/// Since all edge data is one big text block, we need to
/// accumulate all TINYXML_TEXT data and then parse it. This
/// approach effectively skips all comments or other node
......@@ -174,29 +174,29 @@ namespace Nektar
/// missing element numbers due to the text block format.
std::string edgeStr;
int indx;
while(edge)
{
int err = edge->QueryIntAttribute("ID",&indx);
ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
TiXmlNode *child = edge->FirstChild();
edgeStr.clear();
if (child->Type() == TiXmlNode::TINYXML_TEXT)
{
edgeStr += child->ToText()->ValueStr();
}
/// Now parse out the edges, three fields at a time.
int vertex1, vertex2;
std::istringstream edgeDataStrm(edgeStr.c_str());
try
{
while (!edgeDataStrm.fail())
{
edgeDataStrm >> vertex1 >> vertex2;
// Must check after the read because we may be
// at the end and not know it. If we are at
// the end we will add a duplicate of the last
......@@ -204,11 +204,11 @@ namespace Nektar
if (!edgeDataStrm.fail())
{
PointGeomSharedPtr vertices[2] = {GetVertex(vertex1), GetVertex(vertex2)};
SegGeomSharedPtr edge;
it = m_curvedEdges.find(indx);
if(it == m_curvedEdges.end())
{
edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices);
......@@ -219,7 +219,7 @@ namespace Nektar
edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices, it->second);
edge->SetGlobalID(indx); //Set global mesh id
}
m_segGeoms[indx] = edge;
}
}
......@@ -228,7 +228,7 @@ namespace Nektar
{
NEKERROR(ErrorUtil::efatal, (std::string("Unable to read edge data: ") + edgeStr).c_str());
}
edge = edge->NextSiblingElement("E");
}
}
......@@ -257,14 +257,14 @@ namespace Nektar
while (element)
{
std::string elementType(element->ValueStr());
ASSERTL0(elementType == "Q" || elementType == "T",
(std::string("Unknown 2D element type: ") + elementType).c_str());
string IsCompressed;
element->QueryStringAttribute("COMPRESSED",&IsCompressed);
if(IsCompressed.size())
element->QueryStringAttribute("COMPRESSED",&IsCompressed);
if(IsCompressed.size())
{
ASSERTL0(boost::iequals(IsCompressed,
LibUtilities::CompressData:: GetCompressString()),
......@@ -359,7 +359,7 @@ namespace Nektar
};
QuadGeomSharedPtr quadgeom;
if (it == m_curvedFaces.end())
if (it == m_curvedEdges.end())
{
quadgeom =
MemoryManager<QuadGeom>::AllocateSharedPtr(
......@@ -383,9 +383,9 @@ namespace Nektar
int indx;
int err = element->QueryIntAttribute("ID", &indx);
ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
it = m_curvedFaces.find(indx);
/// Read text element description.
TiXmlNode* elementChild = element->FirstChild();
std::string elementStr;
......@@ -787,7 +787,7 @@ namespace Nektar
{
triGeomShPtr = boost::dynamic_pointer_cast<TriGeom>(*geomIter);
quadGeomShPtr = boost::dynamic_pointer_cast<QuadGeom>(*geomIter);
if (triGeomShPtr || quadGeomShPtr)
{
int edgeNum;
......@@ -963,7 +963,7 @@ namespace Nektar
const LibUtilities::PointsKey pkey(numpoints,expansion->m_basisKeyVector[edge_id].GetPointsType());
return LibUtilities::BasisKey(expansion->m_basisKeyVector[edge_id].GetBasisType(),nummodes,pkey);
}
ASSERTL0(false, "Unable to determine edge points type.");
return LibUtilities::NullBasisKey;
}
......
......@@ -514,8 +514,6 @@ namespace Nektar
TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
TiXmlElement* field = NULL;
CurveMap::iterator it;
/// Look for elements in ELEMENT block.
field = mesh->FirstChildElement("ELEMENT");
......@@ -575,11 +573,8 @@ namespace Nektar
boost::static_pointer_cast<TriGeom>(face);
}
it = m_curvedVolumes.find(indx);
TetGeomSharedPtr tetgeom(
MemoryManager<TetGeom>::AllocateSharedPtr(
tfaces, it == m_curvedVolumes.end() ?
CurveSharedPtr() : it->second));
TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>
::AllocateSharedPtr(tfaces));
tetgeom->SetGlobalID(indx);
m_tetGeoms[indx] = tetgeom;
PopulateFaceToElMap(tetgeom, 4);
......@@ -788,11 +783,7 @@ namespace Nektar
ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
it = m_curvedVolumes.find(indx);
TetGeomSharedPtr tetgeom(
MemoryManager<TetGeom>::AllocateSharedPtr(
tfaces, it == m_curvedVolumes.end() ?
CurveSharedPtr() : it->second));
TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>::AllocateSharedPtr(tfaces));
tetgeom->SetGlobalID(indx);
m_tetGeoms[indx] = tetgeom;
......
......@@ -57,12 +57,10 @@ namespace Nektar
m_shapeType = LibUtilities::eTetrahedron;
}
TetGeom::TetGeom(const TriGeomSharedPtr faces[],
CurveSharedPtr curve) :
TetGeom::TetGeom(const TriGeomSharedPtr faces[]) :
Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
{
m_shapeType = LibUtilities::eTetrahedron;
m_curve = curve;
/// Copy the face shared pointers
m_faces.insert(m_faces.begin(), faces, faces+TetGeom::kNfaces);
......
......@@ -44,16 +44,15 @@ namespace Nektar
{
namespace SpatialDomains
{
class TetGeom: public LibUtilities::GraphVertexObject,
class TetGeom: public LibUtilities::GraphVertexObject,
public Geometry3D
{
public:
SPATIAL_DOMAINS_EXPORT TetGeom();
SPATIAL_DOMAINS_EXPORT TetGeom(
const TriGeomSharedPtr faces[],
CurveSharedPtr curve = CurveSharedPtr());
SPATIAL_DOMAINS_EXPORT TetGeom ();
SPATIAL_DOMAINS_EXPORT TetGeom (const TriGeomSharedPtr faces[]);
SPATIAL_DOMAINS_EXPORT ~TetGeom();
SPATIAL_DOMAINS_EXPORT static const int kNverts = 4;
SPATIAL_DOMAINS_EXPORT static const int kNedges = 6;
SPATIAL_DOMAINS_EXPORT static const int kNqfaces = 0;
......
......@@ -637,11 +637,8 @@ void OutputNekpp::WriteXmlCurves(TiXmlElement *pRoot)
}
}
}
if (!curve)
{
return;
}
TiXmlElement *curved = new TiXmlElement("CURVED");
......@@ -695,32 +692,28 @@ void OutputNekpp::WriteXmlCurves(TiXmlElement *pRoot)
m_mesh->m_spaceDim >= 2)
{
vector<ElementSharedPtr>::iterator it;
int d = m_mesh->m_expDim;
std::string tag = d == 2 ? "F" : "V";
std::string attr = d == 2 ? "FACEID" : "VOLID";
for (it = m_mesh->m_element[d].begin();
it != m_mesh->m_element[d].end();
for (it = m_mesh->m_element[m_mesh->m_expDim].begin();
it != m_mesh->m_element[m_mesh->m_expDim].end();
++it)
{
// Only generate face curve if there are volume nodes
if ((*it)->GetVolumeNodes().size() > 0)
{
TiXmlElement *e = new TiXmlElement(tag);
TiXmlElement *e = new TiXmlElement("F");
e->SetAttribute("ID", facecnt++);
e->SetAttribute(attr, (*it)->GetId());
e->SetAttribute("FACEID", (*it)->GetId());
e->SetAttribute("NUMPOINTS", (*it)->GetNodeCount());
e->SetAttribute(
"TYPE",
LibUtilities::kPointsTypeStr[(*it)->GetCurveType()]);
TiXmlText *t0 = new TiXmlText((*it)->GetXmlCurveString());
e->LinkEndChild(t0);
curved->LinkEndChild(e);
}
}
}
if (m_mesh->m_expDim == 3)
else if (m_mesh->m_expDim == 3)
{
FaceSet::iterator it2;
for (it2 = m_mesh->m_faceSet.begin();
......@@ -848,7 +841,7 @@ void OutputNekpp::WriteXmlCurves(TiXmlElement *pRoot)
}
}
// 2D elements in 3-space, output face curvature information
else if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim >= 2)
else if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
{
vector<ElementSharedPtr>::iterator it;
for (it = m_mesh->m_element[m_mesh->m_expDim].begin();
......@@ -1098,7 +1091,7 @@ void OutputNekpp::WriteXmlExpansions(TiXmlElement *pRoot)
if (m_mesh->m_fields.size() == 0)
{
exp->SetAttribute("FIELDS", "u,v,w");
exp->SetAttribute("FIELDS", "u");
}
else
{
......
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