Commit 7e7504e5 authored by Dave Moxey's avatar Dave Moxey

Add changes from Chris

parent 9f319115
......@@ -716,8 +716,6 @@ namespace Nektar
Array<OneD, int> vwgt(nWeight, 1);
Array<OneD, int> vsize(nGraphVerts, 1);
vcnt = 0;
for ( boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
vertit != vertit_end;
++vertit)
......
......@@ -1441,6 +1441,20 @@ namespace Nektar
ASSERTL1(cnt == nbcs,"Failed to visit all boundary condtiions");
}
/**
* @brief Reset this field, so that geometry information can be updated.
*/
void DisContField1D::v_Reset()
{
ExpList::v_Reset();
// Reset boundary condition expansions.
for (int n = 0; n < m_bndCondExpansions.num_elements(); ++n)
{
m_bndCondExpansions[n]->Reset();
}
}
/**
* Search through the edge expansions and identify which ones
* have Robin/Mixed type boundary conditions. If find a Robin
......
......@@ -240,7 +240,8 @@ namespace Nektar
virtual void v_GetBoundaryToElmtMap(
Array<OneD,int> &ElmtID, Array<OneD,int> &VertID);
virtual void v_Reset();
/// Evaluate all boundary conditions at a given time..
virtual void v_EvaluateBoundaryConditions(
const NekDouble time = 0.0,
......
......@@ -536,50 +536,49 @@ namespace Nektar
int i,j,k;
int nEdgeCoeffs;
for(int i = 0; i < m_coordim; ++i)
if (m_curve)
{
if (!m_curve)
{
continue;
}
int npts = m_curve->m_points.size();
int nEdgePts = (int)sqrt(static_cast<NekDouble>(npts));
Array<OneD,NekDouble> tmp(npts);
LibUtilities::PointsKey curveKey(nEdgePts, m_curve->m_ptype);
Array<OneD, NekDouble> tmp(npts);
Array<OneD, NekDouble> tmp2(m_xmap->GetTotPoints());
LibUtilities::PointsKey curveKey(
nEdgePts, m_curve->m_ptype);
// Sanity checks:
// - Curved faces should have square number of points;
// - Each edge should have sqrt(npts) points.
ASSERTL0(nEdgePts*nEdgePts == npts,
ASSERTL0(nEdgePts * nEdgePts == npts,
"NUMPOINTS should be a square number in"
" quadrilteral "
+ boost::lexical_cast<string>(m_globalID));
for (j = 0; j < kNedges; ++j)
for (i = 0; i < kNedges; ++i)
{
ASSERTL0(
m_edges[j]->GetXmap()->GetNcoeffs() == nEdgePts,
m_edges[i]->GetXmap()->GetNcoeffs() == nEdgePts,
"Number of edge points does not correspond to "
"number of face points in quadrilateral "
+ boost::lexical_cast<string>(m_globalID));
}
for (j = 0; j < npts; ++j)
for (i = 0; i < m_coordim; ++i)
{
tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
}
for (j = 0; j < npts; ++j)
{
tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
}
// Interpolate m_curve points to GLL points
Array<OneD, NekDouble> tmp2(m_xmap->GetTotPoints());
LibUtilities::Interp2D(
curveKey, curveKey, tmp,
m_xmap->GetBasis(0)->GetPointsKey(),
m_xmap->GetBasis(1)->GetPointsKey(),
tmp2);
// Forwards transform to get coefficient space.
m_xmap->FwdTrans(tmp2, m_coeffs[i]);
// Interpolate m_curve points to GLL points
LibUtilities::Interp2D(
curveKey, curveKey, tmp,
m_xmap->GetBasis(0)->GetPointsKey(),
m_xmap->GetBasis(1)->GetPointsKey(),
tmp2);
// Forwards transform to get coefficient space.
m_xmap->FwdTrans(tmp2, m_coeffs[i]);
}
}
// Now fill in edges.
......
......@@ -370,6 +370,11 @@ namespace Nektar
NEKERROR(ErrorUtil::ewarning, err.c_str());
}
LibUtilities::PointsKey fkey(npts,m_curve->m_ptype);
DNekMatSharedPtr I0 =
LibUtilities::PointsManager()[fkey]->GetI(pkey);
NekVector<NekDouble> out(npts+1);
for(int i = 0; i < m_coordim; ++i)
{
// Load up coordinate values into tmp
......@@ -379,13 +384,7 @@ namespace Nektar
}
// Interpolate to GLL points
DNekMatSharedPtr I0;
LibUtilities::PointsKey fkey(npts,m_curve->m_ptype);
I0 = LibUtilities::PointsManager()[fkey]->GetI(pkey);
NekVector<NekDouble> in (npts, tmp, eWrapper);
NekVector<NekDouble> out(npts+1);
out = (*I0)*in;
m_xmap->FwdTrans(out.GetPtr(), m_coeffs[i]);
......
......@@ -724,7 +724,7 @@ namespace Nektar
SetUpXmap();
SetUpCoeffs(m_xmap->GetNcoeffs());
}
/**
* @brief Set up the #m_xmap object by determining the order of each
* direction from derived faces.
......@@ -734,17 +734,17 @@ namespace Nektar
vector<int> tmp;
tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
int order0 = *max_element(tmp.begin(), tmp.end());
tmp.clear();
tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(0));
int points0 = *max_element(tmp.begin(), tmp.end());
tmp.clear();
tmp.push_back(order0);
tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
int order1 = *max_element(tmp.begin(), tmp.end());
tmp.clear();
tmp.push_back(points0);
tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(1));
......
......@@ -456,19 +456,15 @@ namespace Nektar
int i,j,k;
int nEdgeCoeffs = m_xmap->GetEdgeNcoeffs(0);
for(int i = 0; i < m_coordim; ++i)
if (m_curve)
{
if (!m_curve)
{
continue;
}
int pdim = LibUtilities::PointsManager()[
LibUtilities::PointsKey(2, m_curve->m_ptype)]
->GetPointsDim();
// Deal with 2D points type separately (e.g. electrostatic
// or Fekete points) to 1D tensor product.
// Deal with 2D points type separately
// (e.g. electrostatic or Fekete points) to 1D tensor
// product.
if (pdim == 2)
{
int N = m_curve->m_points.size();
......@@ -480,16 +476,15 @@ namespace Nektar
" triangle "
+ boost::lexical_cast<string>(m_globalID));
for (int j = 0; j < kNedges; ++j)
for (i = 0; i < kNedges; ++i)
{
ASSERTL0(
m_edges[j]->GetXmap()->GetNcoeffs() == nEdgePts,
"Number of edge points does not correspond to "
"number of face points in triangle "
m_edges[i]->GetXmap()->GetNcoeffs() == nEdgePts,
"Number of edge points does not correspond "
"to number of face points in triangle "
+ boost::lexical_cast<string>(m_globalID));
}
// Create a StdNodalTriExp.
const LibUtilities::PointsKey P0(
nEdgePts, LibUtilities::eGaussLobattoLegendre);
const LibUtilities::PointsKey P1(
......@@ -498,37 +493,41 @@ namespace Nektar
LibUtilities::eOrtho_A, nEdgePts, P0);
const LibUtilities::BasisKey T1(
LibUtilities::eOrtho_B, nEdgePts, P1);
StdRegions::StdNodalTriExpSharedPtr t =
MemoryManager<StdRegions::StdNodalTriExp>
::AllocateSharedPtr(T0, T1, m_curve->m_ptype);
Array<OneD, NekDouble> phys(
max(t->GetTotPoints(), m_xmap->GetTotPoints()));
max(nEdgePts*nEdgePts, m_xmap->GetTotPoints()));
Array<OneD, NekDouble> tmp(nEdgePts*nEdgePts);
for (int j = 0; j < N; ++j)
for(i = 0; i < m_coordim; ++i)
{
phys[j] = (m_curve->m_points[j]->GetPtr())[i];
// Create a StdNodalTriExp.
StdRegions::StdNodalTriExpSharedPtr t =
MemoryManager<StdRegions::StdNodalTriExp>
::AllocateSharedPtr(T0, T1, 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::Interp2D(
P0, P1, tmp,
m_xmap->GetBasis(0)->GetPointsKey(),
m_xmap->GetBasis(1)->GetPointsKey(),
phys);
// Forwards transform to get coefficient space.
m_xmap->FwdTrans(phys, m_coeffs[i]);
}
Array<OneD, NekDouble> tmp(nEdgePts*nEdgePts);
t->BwdTrans(phys, tmp);
// Interpolate points to standard region.
LibUtilities::Interp2D(
P0, P1, tmp,
m_xmap->GetBasis(0)->GetPointsKey(),
m_xmap->GetBasis(1)->GetPointsKey(),
phys);
// Forwards transform to get coefficient space.
m_xmap->FwdTrans(phys, m_coeffs[i]);
}
else if (pdim == 1)
{
int npts = m_curve->m_points.size();
int nEdgePts = (int)sqrt(static_cast<NekDouble>(npts));
Array<OneD,NekDouble> tmp(npts);
Array<OneD, NekDouble> tmp (npts);
Array<OneD, NekDouble> phys(m_xmap->GetTotPoints());
LibUtilities::PointsKey curveKey(
nEdgePts, m_curve->m_ptype);
......@@ -540,30 +539,33 @@ namespace Nektar
" triangle "
+ boost::lexical_cast<string>(m_globalID));
for (int j = 0; j < kNedges; ++j)
for (i = 0; i < kNedges; ++i)
{
ASSERTL0(
m_edges[j]->GetXmap()->GetNcoeffs() == nEdgePts,
m_edges[i]->GetXmap()->GetNcoeffs() == nEdgePts,
"Number of edge points does not correspond to "
"number of face points in triangle "
+ boost::lexical_cast<string>(m_globalID));
}
for (int j = 0; j < npts; ++j)
for (i = 0; i < m_coordim; ++i)
{
tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
for (j = 0; j < npts; ++j)
{
tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
}
// Interpolate curve points to standard triangle
// points.
LibUtilities::Interp2D(
curveKey, curveKey, tmp,
m_xmap->GetBasis(0)->GetPointsKey(),
m_xmap->GetBasis(1)->GetPointsKey(),
phys);
// Forwards transform to get coefficient space.
m_xmap->FwdTrans(phys, m_coeffs[i]);
}
// Interpolate curve points to standard triangle points.
Array<OneD, NekDouble> phys(m_xmap->GetTotPoints());
LibUtilities::Interp2D(
curveKey, curveKey, tmp,
m_xmap->GetBasis(0)->GetPointsKey(),
m_xmap->GetBasis(1)->GetPointsKey(),
phys);
// Forwards transform to get coefficient space.
m_xmap->FwdTrans(phys, m_coeffs[i]);
}
else
{
......
......@@ -42,18 +42,24 @@
namespace Nektar
{
class CoupledAssemblyMap : public MultiRegions::AssemblyMapCG
{
public:
CoupledAssemblyMap(
const LibUtilities::SessionReaderSharedPtr &pSession,
const SpatialDomains::MeshGraphSharedPtr &graph,
const MultiRegions::AssemblyMapCGSharedPtr &cgMap,
const SpatialDomains::BoundaryConditionsSharedPtr &boundConds,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields);
};
typedef boost::shared_ptr<CoupledAssemblyMap> CoupledAssemblyMapSharedPtr;
/**
* @brief Modified version of MultiRegions::AssemblyMapCG that allows for
* coupled fields [u,v,w] instead of individual scalar fields u, v and w.
*/
class CoupledAssemblyMap : public MultiRegions::AssemblyMapCG
{
public:
CoupledAssemblyMap(
const LibUtilities::SessionReaderSharedPtr &pSession,
const SpatialDomains::MeshGraphSharedPtr &graph,
const MultiRegions::AssemblyMapCGSharedPtr &cgMap,
const SpatialDomains::BoundaryConditionsSharedPtr &boundConds,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields);
};
typedef boost::shared_ptr<CoupledAssemblyMap> CoupledAssemblyMapSharedPtr;
}
#endif
......@@ -43,44 +43,50 @@ using namespace Nektar::SolverUtils;
namespace Nektar
{
class IterativeElasticSystem : public LinearElasticSystem
/**
* @brief Class for iterative elastic system, in which linear elasticity is
* applied in substeps to attain a large deformation.
*/
class IterativeElasticSystem : public LinearElasticSystem
{
public:
/// Class may only be instantiated through the MemoryManager.
friend class MemoryManager<IterativeElasticSystem>;
/// Creates an instance of this class
static EquationSystemSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession)
{
public:
/// Class may only be instantiated through the MemoryManager.
friend class MemoryManager<IterativeElasticSystem>;
EquationSystemSharedPtr p = MemoryManager<
IterativeElasticSystem>::AllocateSharedPtr(pSession);
p->InitObject();
return p;
}
/// Creates an instance of this class
static EquationSystemSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession)
{
EquationSystemSharedPtr p = MemoryManager<
IterativeElasticSystem>::AllocateSharedPtr(pSession);
p->InitObject();
return p;
}
/// Name of class
static std::string className;
/// Name of class
static std::string className;
protected:
/// Number of steps to split deformation across.
int m_numSteps;
/// Flag determining whether to repeat boundary conditions.
bool m_repeatBCs;
/// Storage for boundary conditions.
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_savedBCs;
/// Vector of boundary regions to deform.
vector<int> m_toDeform;
protected:
/// Number of steps to split deformation across.
int m_numSteps;
/// Flag determining whether to repeat boundary conditions.
bool m_repeatBCs;
/// Storage for boundary conditions.
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_savedBCs;
/// Vector of boundary regions to deform.
vector<int> m_toDeform;
IterativeElasticSystem(
const LibUtilities::SessionReaderSharedPtr& pSession);
IterativeElasticSystem(
const LibUtilities::SessionReaderSharedPtr& pSession);
virtual void v_InitObject();
virtual void v_GenerateSummary(SolverUtils::SummaryList& s);
virtual void v_DoSolve();
virtual void v_InitObject();
virtual void v_GenerateSummary(SolverUtils::SummaryList& s);
virtual void v_DoSolve();
void WriteGeometry(const int i);
};
void WriteGeometry(const int i);
};
}
#endif
......@@ -43,72 +43,79 @@ using namespace Nektar::SolverUtils;
namespace Nektar
{
class LinearElasticSystem : public EquationSystem
/**
* @brief Base class for linear elastic system.
*/
class LinearElasticSystem : public EquationSystem
{
public:
/// Class may only be instantiated through the MemoryManager.
friend class MemoryManager<LinearElasticSystem>;
/// Creates an instance of this class
static EquationSystemSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession)
{
public:
/// Class may only be instantiated through the MemoryManager.
friend class MemoryManager<LinearElasticSystem>;
EquationSystemSharedPtr p = MemoryManager<
LinearElasticSystem>::AllocateSharedPtr(pSession);
p->InitObject();
return p;
}
/// Creates an instance of this class
static EquationSystemSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession)
{
EquationSystemSharedPtr p = MemoryManager<
LinearElasticSystem>::AllocateSharedPtr(pSession);
p->InitObject();
return p;
}
/// Name of class
static std::string className;
/// Name of class
static std::string className;
void BuildMatrixSystem();
void BuildMatrixSystem();
void SetStaticCondBlock(
const int n,
const LocalRegions::ExpansionSharedPtr exp,
Array<TwoD, DNekMatSharedPtr> &mat);
void SetStaticCondBlock(
const int n,
const LocalRegions::ExpansionSharedPtr exp,
Array<TwoD, DNekMatSharedPtr> &mat);
DNekMatSharedPtr BuildLaplacianIJMatrix(
const int k1,
const int k2,
const NekDouble scale,
LocalRegions::ExpansionSharedPtr exp);
DNekMatSharedPtr BuildLaplacianIJMatrix(
const int k1,
const int k2,
const NekDouble scale,
LocalRegions::ExpansionSharedPtr exp);
protected:
LinearElasticSystem(
const LibUtilities::SessionReaderSharedPtr& pSession);
protected:
LinearElasticSystem(
const LibUtilities::SessionReaderSharedPtr& pSession);
/// Poisson ratio.
NekDouble m_nu;
/// Young's modulus.
NekDouble m_E;
/// Parameter dictating amount of thermal stress to add.
NekDouble m_beta;
/// Assembly map for the coupled (u,v,w) system.
CoupledAssemblyMapSharedPtr m_assemblyMap;
/// Schur complement boundary-boundary matrix.
DNekScalBlkMatSharedPtr m_schurCompl;
/// Matrix of elemental \f$ B^{-1}D \f$ components
DNekScalBlkMatSharedPtr m_BinvD;
/// Matrix of elemental \f$ C \f$ components
DNekScalBlkMatSharedPtr m_C;
/// Matrix of elemental \f$ D^{-1} \f$ components
DNekScalBlkMatSharedPtr m_Dinv;
/// Boundary maps for each of the fields.
Array<OneD, Array<OneD, unsigned int> > m_bmap;
/// Interior maps for each of the fields.
Array<OneD, Array<OneD, unsigned int> > m_imap;
/// Storage for the temperature terms.
Array<OneD, Array<OneD, NekDouble> > m_temperature;
/// Storage for the thermal stress terms.
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_stress;
/// Poisson ratio.
NekDouble m_nu;
/// Young's modulus
NekDouble m_E;
/// Parameter dictating amount of thermal stress to add.
NekDouble m_beta;
/// Assembly map for the coupled (u,v,w) system.
CoupledAssemblyMapSharedPtr m_assemblyMap;
/// Schur complement boundary-boundary matrix.
DNekScalBlkMatSharedPtr m_schurCompl;
/// Matrix of elemental \f$ B^{-1}D \f$ components
DNekScalBlkMatSharedPtr m_BinvD;
/// Matrix of elemental \f$ C \f$ components
DNekScalBlkMatSharedPtr m_C;
/// Matrix of elemental \f$ D^{-1} \f$ components
DNekScalBlkMatSharedPtr m_Dinv;
Array<OneD, Array<OneD, unsigned int> > m_bmap;
Array<OneD, Array<OneD, unsigned int> > m_imap;
/// Storage for the temperature terms.
Array<OneD, Array<OneD, NekDouble> > m_temperature;
/// Storage for the thermal stress terms.
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_stress;
virtual void v_InitObject();
virtual void v_GenerateSummary(SolverUtils::SummaryList& s);
virtual void v_DoSolve();
virtual void v_ExtraFldOutput(
std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
std::vector<std::string> &variables);
};
virtual void v_InitObject();
virtual void v_GenerateSummary(SolverUtils::SummaryList& s);
virtual void v_DoSolve();
virtual void v_ExtraFldOutput(
std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
std::vector<std::string> &variables);
};
}
#endif
......@@ -813,11 +813,6 @@ namespace Nektar
}
LibUtilities::ShapeType surfElType = surfEl->GetConf().m_e;
// Hack
if (surfElType == LibUtilities::eQuadrilateral)
{
surfElType = LibUtilities::eTriangle;
}
if (!found)
{
......
......@@ -402,7 +402,8 @@ namespace Nektar
*
* The last step is to eliminate duplicate edges/faces and reenumerate.
*
* NOTE: This routine does not copy high-order information yet!
* NOTE: This routine does not copy face-interior high-order information