Commit d3bfd923 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/FieldConvertEquiOutput' into 'master'

Feature/field convert equi output

See merge request !421
parents 7ab1dd08 6e9bfedc
......@@ -97,11 +97,13 @@ Specifically, FieldConvert has these additional functionalities
\item \inltt{C0Projection}: Computes the C0 projection of a given output file;
\item \inltt{QCriterion}: Computes the Q-Criterion for a given output file;
\item \inltt{concatenate}: Concatenate a \nekpp binary output (.chk or .fld) field file into single file;
\item \inltt{equispacedoutput}: Write data as equi-spaced output using simplices to represent the data for connecting points;
\item \inltt{extract}: Extract a boundary field;
\item \inltt{interpfield}: Interpolates one field to another, requires fromxml, fromfld to be defined;
\item \inltt{interppointdatatofld}: Interpolates given discrete data using a finite difference approximation to a fld file given an xml file;
\item \inltt{interppoints}: Interpolates a set of points to another, requires fromfld and fromxml to be defined, a line or plane of points can be defined;
\item \inltt{scaleinputfld}: rescale input field by a constant factor.
\item \inltt{isocontour}: Extract an isocontour of ``fieldid'' variable and at value ``fieldvalue''. Optionally ``fieldstr'' can be specified for a string defiition or ``smooth'' for smoothing;
\item \inltt{scaleinputfld}: Rescale input field by a constant factor.
\item \inltt{vorticity}: Computes the vorticity field.
\end{enumerate}
The module list above can be seen by running the command
......@@ -174,6 +176,24 @@ to visualise either in Tecplot or in Paraview the result.
%
%
%
\subsubsection{Equi-spaced output of data: \textit{equispacedoutput} module}
This module interpolates the output data to an truly equispaced set of
points (not equispaced along the collapsed coordinate
system). Therefore a tetrahedron is represented by a tetrahedral
number of poinst. This produces much smaller output files. The points
are then connected together by simplices (triangles and tetrahedrons).
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m equispacedouptput test.xml test.fld test-boundary.dat
\end{lstlisting}
\begin{notebox}
Currently this option is only set up for triangles, quadrilaterals,
tetrahedrons and prisms. It also only is currently used in tecplot
output.
\end{notebox}
\subsubsection{Extract a boundary region: \textit{extract} module}
The boundary region of a domain can be extracted from the output
data using the following command line
......@@ -333,12 +353,49 @@ where \inltt{npts1,npts2} is the number of equispaced points in each
direction and $(x0,y0,z0)$, $(x1,y1,z1)$, $(x2,y2,z2)$ and $(x3,y3,z3)$
define the plane of points specified in a clockwise or anticlockwise direction.
%
\begin{notebox}
\begin{notebox}
This module does not run in parallel.
\end{notebox}
%
%
%
\subsubsection{Isoncontour extraction: \textit{iscontour} module}
Extract an isocontour from a field file. This option automatically
take the field to an equispaced distribution of points connected by
linear simplicies of triangles or tetrahedrons. The linear simplices
are then inspected to extract the isocontour of interest. To specify
the field \inltt{fieldid} can be provided giving the id of the field
of interest and \inltt{fieldvalue} provides the value of the
isocontour to be extracted.
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m isocontour:fieldid=2:fieldvalue=0.5 test.xml test.fld \
test-isocontour.dat
\end{lstlisting}
Alternatively \inltt{fieldstr}=''u+v'' can be specified to calculate
the field $u^2$ and extract its isocontour. You can also specify
\inltt{fieldname}=''UplusV'' to define the name of the isocontour in
the .dat file, i.e.
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m isocontour:fieldstr="u+v":fieldvalue=0.5:\
fieldname="UplusV" test.xml test.fld test-isocontour.dat
\end{lstlisting}
Optionally \inltt{smooth} can be specified to smooth the isocontour
with default values \inltt{smoothnegdiffusion}=0.495,
\inltt{smoothnegdiffusion}=0.5 and \inltt{smoothiter}=100. This option
automatically calls a \inltt{globalcondense} option which remove
multiply defined verties from the simplex definition which arise as
isocontour are generated element by element.
\begin{notebox}
Currently this option is only set up for triangles, quadrilaterals,
tetrahedrons and prisms.
\end{notebox}
\subsubsection{Scale a given .fld: \textit{scaleinputfld} module}
\begin{lstlisting}[style=BashInputStyle]
......
......@@ -342,6 +342,41 @@ namespace Nektar
return returnval;
}
inline int GetNumberOfCoefficients(ShapeType shape, int na, int nb, int nc)
{
int returnval = 0;
switch(shape)
{
case eSegment:
returnval = na;
break;
case eTriangle:
returnval = StdTriData::getNumberOfCoefficients(na,nb);
break;
case eQuadrilateral:
returnval = na*nb;
break;
case eTetrahedron:
returnval = StdTetData::getNumberOfCoefficients(na,nb,nc);
break;
case ePyramid:
returnval = StdPyrData::getNumberOfCoefficients(na,nb,nc);
break;
case ePrism:
returnval = StdPrismData::getNumberOfCoefficients(na,nb,nc);
break;
case eHexahedron:
returnval = na*nb*nc;
break;
default:
ASSERTL0(false,"Unknown Shape Type");
break;
}
return returnval;
}
}
}
......
......@@ -83,7 +83,7 @@ namespace Nektar
func(PFD3 p) : func3(p), size(3) {};
func(PFD4 p) : func4(p), size(4) {};
union // Pointer to a function
union // Pointer to a function
{
PFD1 func1;
PFD2 func2;
......@@ -120,7 +120,7 @@ namespace Nektar
{
functions()
{
// Add all of the functions from math.h
// Add all of the functions from math.h
add
("abs", std::abs)
("asin", asin)
......@@ -191,7 +191,7 @@ namespace Nektar
| (root_node_d[ch_p('-')] >> factor);
parameter = leaf_node_d[ lexeme_d[
(alpha_p | '_' | '$') >> *(alnum_p | '_' | '$')
(alpha_p | '_' | '$') >> *(alnum_p | '_' | '$')
] ] >> op;
function = root_node_d[functions_p] >>
......@@ -335,7 +335,7 @@ namespace Nektar
if(m_constant[it->second] != value)
{
std::string errormsg("Attempt to add numerically different constants under the same name: ");
errormsg += name;
errormsg += name;
std::cout << errormsg << std::endl;
}
//ASSERTL1(m_constant[it->second] == value, "Attempt to add numerically different constants under the same name: " + name);
......@@ -581,6 +581,7 @@ namespace Nektar
}
// wrapper function call
void AnalyticExpressionEvaluator::Evaluate(
const int expression_id,
const Array<OneD, const NekDouble>& x,
......@@ -591,9 +592,34 @@ namespace Nektar
{
m_timer.Start();
const int num_points = x.num_elements();
ASSERTL1(m_executionStack.size() > expression_id, "unknown analytic expression, it must first be defined with DefineFunction(...)");
ASSERTL1(result.num_elements() >= num_points, "destination array must have enough capacity to store expression values at each given point");
std::vector<Array<OneD, const NekDouble> >points;
points.push_back(x);
points.push_back(y);
points.push_back(z);
points.push_back(t);
Evaluate(expression_id,points,result);
m_timer.Stop();
m_total_eval_time += m_timer.TimePerTest(1);
}
void AnalyticExpressionEvaluator::Evaluate(
const int expression_id,
const std::vector<Array<OneD, const NekDouble> > points,
Array<OneD, NekDouble>& result)
{
m_timer.Start();
const int num_points = points[0].num_elements();
ASSERTL1(m_executionStack.size() > expression_id,
"unknown analytic expression, it must first be defined "
"with DefineFunction(...)");
ASSERTL1(result.num_elements() >= num_points,
"destination array must have enough capacity to store "
"expression values at each given point");
ExecutionStack &stack = m_executionStack[expression_id];
......@@ -602,16 +628,16 @@ namespace Nektar
/// Lets split the work into cache-sized chunks.
/// Ahtung, magic constant!
const int max_chunk_size = 1024;
/// please don't remove brackets around std::min, it screws up windows compilation
const int nvals = points.size();
const int chunk_size = (std::min)(max_chunk_size, num_points);
if (m_state.size() < chunk_size * m_state_sizes[expression_id] )
if (m_state.size() < chunk_size * m_state_sizes[expression_id])
{
m_state.resize( m_state_sizes[expression_id] * chunk_size, 0.0 );
m_state.resize(m_state_sizes[expression_id] * chunk_size, 0.0);
}
if (m_variable.size() < 4 * chunk_size )
if (m_variable.size() < nvals * chunk_size)
{
m_variable.resize( 4 * chunk_size, 0.0);
m_variable.resize( nvals * chunk_size, 0.0);
}
if (result.num_elements() < num_points)
{
......@@ -625,10 +651,10 @@ namespace Nektar
const int this_chunk_size = (std::min)(work_left, 1024);
for (int i = 0; i < this_chunk_size; i++)
{
m_variable[i+this_chunk_size*0] = x[offset + i];
m_variable[i+this_chunk_size*1] = y[offset + i];
m_variable[i+this_chunk_size*2] = z[offset + i];
m_variable[i+this_chunk_size*3] = t[offset + i];
for(int j = 0; j < nvals; ++j)
{
m_variable[i+this_chunk_size*j] = points[j][offset + i];
}
}
for (int i = 0; i < stack.size(); i++)
{
......@@ -646,48 +672,6 @@ namespace Nektar
}
void AnalyticExpressionEvaluator::EvaluateAtPoints(
const int expression_id,
const std::vector<Array<OneD, const NekDouble> > points,
Array<OneD, NekDouble>& result)
{
m_timer.Start();
/// \todo test this function properly/update as the method above
ASSERTL1(m_executionStack.size() > expression_id, "unknown analytic expression, it must first be defined with DefineFunction(...)");
ExecutionStack& stack = m_executionStack[expression_id];
VariableMap& variableMap = m_stackVariableMap[expression_id];
const int num = points[0].num_elements();
m_state.resize(m_state_sizes[expression_id]*num);
// assuming all points have same # of coordinates
m_variable.resize(4*num,0.0);
for (int i = 0; i < points.size(); i++)
{
for (VariableMap::const_iterator it = variableMap.begin(); it != variableMap.end(); ++it)
{
m_variable[it->second] = points[i][it->second];
}
}
for (int j = 0; j < stack.size(); j++)
{
(*stack[j]).run_many(num);
}
for (int i = 0; i < num; ++i)
{
result[i] = m_state[i];
}
m_timer.Stop();
m_total_eval_time += m_timer.TimePerTest(1);
}
AnalyticExpressionEvaluator::PrecomputedValue AnalyticExpressionEvaluator::PrepareExecutionAsYouParse(
const ParsedTreeIterator& location,
ExecutionStack& stack,
......@@ -805,7 +789,7 @@ namespace Nektar
case E_COSH:
stack.push_back ( makeStep<EvalCosh>( stateIndex, stateIndex ) );
return std::make_pair(false,0);
case E_EXP:
case E_EXP:
stack.push_back ( makeStep<EvalExp>( stateIndex, stateIndex ) );
return std::make_pair(false,0);
case E_FABS:
......@@ -814,7 +798,7 @@ namespace Nektar
case E_FLOOR:
stack.push_back ( makeStep<EvalFloor>( stateIndex, stateIndex ) );
return std::make_pair(false,0);
case E_LOG:
case E_LOG:
stack.push_back ( makeStep<EvalLog>( stateIndex, stateIndex ) );
return std::make_pair(false,0);
case E_LOG10:
......@@ -823,7 +807,7 @@ namespace Nektar
case E_RAD:
stack.push_back ( makeStep<EvalRad>( stateIndex, stateIndex, stateIndex+1 ) );
return std::make_pair(false,0);
case E_SIN:
case E_SIN:
stack.push_back ( makeStep<EvalSin>( stateIndex, stateIndex ) );
return std::make_pair(false,0);
case E_SINH:
......@@ -832,7 +816,7 @@ namespace Nektar
case E_SQRT:
stack.push_back ( makeStep<EvalSqrt>( stateIndex, stateIndex ) );
return std::make_pair(false,0);
case E_TAN:
case E_TAN:
stack.push_back ( makeStep<EvalTan>( stateIndex, stateIndex ) );
return std::make_pair(false,0);
case E_TANH:
......
......@@ -247,9 +247,11 @@ namespace Nektar
const Array<OneD, const NekDouble>&,
Array<OneD, NekDouble>& result);
/// Vectorized evaluation method for expressions depending on unspecified
/// number of variables.
LIB_UTILITIES_EXPORT void EvaluateAtPoints(
LIB_UTILITIES_EXPORT void Evaluate(
const int expression_id,
const std::vector<Array<OneD, const NekDouble> > points,
Array<OneD, NekDouble>& result);
......
This diff is collapsed.
......@@ -45,19 +45,19 @@
namespace Nektar
{
namespace LocalRegions
namespace LocalRegions
{
class PrismExp : virtual public StdRegions::StdPrismExp,
class PrismExp : virtual public StdRegions::StdPrismExp,
virtual public Expansion3D
{
public:
/// \brief Constructor using BasisKey class for quadrature
/// points and order definition
/// points and order definition
LOCAL_REGIONS_EXPORT PrismExp(const LibUtilities::BasisKey &Ba,
const LibUtilities::BasisKey &Bb,
const LibUtilities::BasisKey &Bc,
const SpatialDomains::PrismGeomSharedPtr &geom);
LOCAL_REGIONS_EXPORT PrismExp(const PrismExp &T);
LOCAL_REGIONS_EXPORT ~PrismExp();
......@@ -67,14 +67,14 @@ namespace Nektar
// Integration Methods
//-------------------------------
LOCAL_REGIONS_EXPORT virtual NekDouble v_Integral(
const Array<OneD, const NekDouble>& inarray);
const Array<OneD, const NekDouble>& inarray);
//----------------------------
// Differentiation Methods
//----------------------------
LOCAL_REGIONS_EXPORT virtual void v_PhysDeriv(
const Array<OneD, const NekDouble>& inarray,
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& out_d0,
Array<OneD, NekDouble>& out_d1,
Array<OneD, NekDouble>& out_d2);
......@@ -111,7 +111,7 @@ namespace Nektar
// Evaluation functions
//---------------------------------------
LOCAL_REGIONS_EXPORT virtual void v_GetCoord(
const Array<OneD, const NekDouble> &Lcoords,
const Array<OneD, const NekDouble> &Lcoords,
Array<OneD, NekDouble> &coords);
LOCAL_REGIONS_EXPORT virtual void v_GetCoords(
......@@ -126,7 +126,7 @@ namespace Nektar
LOCAL_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(
const Array<OneD, const NekDouble>& coord,
const Array<OneD, const NekDouble>& physvals);
//---------------------------------------
// Helper functions
......@@ -134,8 +134,8 @@ namespace Nektar
LOCAL_REGIONS_EXPORT virtual int v_GetCoordim();
LOCAL_REGIONS_EXPORT virtual void v_ExtractDataToCoeffs(
const NekDouble *data,
const std::vector<unsigned int > &nummodes,
const int mode_offset,
const std::vector<unsigned int > &nummodes,
const int mode_offset,
NekDouble * coeffs);
LOCAL_REGIONS_EXPORT virtual void v_GetFacePhysVals(
const int face,
......@@ -152,24 +152,24 @@ namespace Nektar
StdRegions::Orientation orient);
LOCAL_REGIONS_EXPORT void v_ComputeFaceNormal(const int face);
//---------------------------------------
// Operator creation functions
//---------------------------------------
LOCAL_REGIONS_EXPORT virtual void v_MassMatrixOp(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(
const int k1,
const int k2,
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
const int k1,
const int k2,
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
......@@ -178,11 +178,6 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_ReduceOrderCoeffs(
int numMin,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
//---------------------------------------
// Matrix creation functions
//---------------------------------------
......@@ -201,10 +196,17 @@ namespace Nektar
LOCAL_REGIONS_EXPORT DNekScalBlkMatSharedPtr CreateStaticCondMatrix(
const MatrixKey &mkey);
LOCAL_REGIONS_EXPORT
virtual void v_GetSimplexEquiSpacedConnectivity(
Array<OneD, int> &conn,
bool standard = true);
private:
LibUtilities::NekManager<MatrixKey, DNekScalMat, MatrixKey::opLess> m_matrixManager;
LibUtilities::NekManager<MatrixKey, DNekScalBlkMat, MatrixKey::opLess> m_staticCondMatrixManager;
LibUtilities::NekManager<MatrixKey, DNekScalMat,
MatrixKey::opLess> m_matrixManager;
LibUtilities::NekManager<MatrixKey, DNekScalBlkMat,
MatrixKey::opLess> m_staticCondMatrixManager;
virtual void v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
......
......@@ -1485,79 +1485,6 @@ namespace Nektar
}
}
void TetExp::v_ReduceOrderCoeffs(
int numMin,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
int n_coeffs = inarray.num_elements();
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nqtot = nquad0*nquad1*nquad2;
int nmodes0 = m_base[0]->GetNumModes();
int nmodes1 = m_base[1]->GetNumModes();
int nmodes2 = m_base[2]->GetNumModes();
int numMax = nmodes0;
Array<OneD, NekDouble> coeff (n_coeffs);
Array<OneD, NekDouble> coeff_tmp1(n_coeffs, 0.0);
Array<OneD, NekDouble> coeff_tmp2(n_coeffs, 0.0);
Array<OneD, NekDouble> phys_tmp(nqtot,0.0);
Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp2,1);
const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
LibUtilities::BasisKey b0 = m_base[0]->GetBasisKey();
LibUtilities::BasisKey b1 = m_base[1]->GetBasisKey();
LibUtilities::BasisKey b2 = m_base[2]->GetBasisKey();
LibUtilities::BasisKey bortho0(
LibUtilities::eOrtho_A, nmodes0, Pkey0);
LibUtilities::BasisKey bortho1(
LibUtilities::eOrtho_B, nmodes1, Pkey1);
LibUtilities::BasisKey bortho2(
LibUtilities::eOrtho_C, nmodes2, Pkey2);
Vmath::Zero(n_coeffs, coeff_tmp2, 1);
StdRegions::StdTetExpSharedPtr m_OrthoTetExp;
StdRegions::StdTetExpSharedPtr m_TetExp;
m_TetExp = MemoryManager<StdRegions::StdTetExp>
::AllocateSharedPtr(b0, b1, b2);
m_OrthoTetExp = MemoryManager<StdRegions::StdTetExp>
::AllocateSharedPtr(bortho0, bortho1, bortho2);
m_TetExp ->BwdTrans(inarray,phys_tmp);
m_OrthoTetExp->FwdTrans(phys_tmp, coeff);
Vmath::Zero(m_ncoeffs,outarray,1);
// filtering
int cnt = 0;
for (int u = 0; u < numMax; ++u)
{
for (int i = 0; i < numMax-u; ++i)
{
Vmath::Vcopy(numMin-u-i,
tmp = coeff+cnt,1,
tmp2 = coeff_tmp1+cnt,1);
cnt += numMax - u - i;
}
}
m_OrthoTetExp->BwdTrans(coeff_tmp1,phys_tmp);
m_TetExp ->FwdTrans(phys_tmp, outarray);
}
void TetExp::v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
......
......@@ -159,11 +159,6 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_ReduceOrderCoeffs(
int numMin,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
LOCAL_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
......
This diff is collapsed.
<
......@@ -134,20 +134,23 @@ namespace Nektar
typedef std::vector<ElementFaceSharedPtr> ElementFaceVector;
typedef boost::shared_ptr<ElementFaceVector> ElementFaceVectorSharedPtr;