本文整理汇总了C++中intrepid::FieldContainer类的典型用法代码示例。如果您正苦于以下问题:C++ FieldContainer类的具体用法?C++ FieldContainer怎么用?C++ FieldContainer使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了FieldContainer类的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: values
void UnitNormalFunction::values(Intrepid::FieldContainer<double> &values, BasisCachePtr basisCache)
{
this->CHECK_VALUES_RANK(values);
int numCells = values.dimension(0);
int numPoints = values.dimension(1);
int spaceDim = basisCache->getSpaceDim();
if (_comp == -1)
{
// check the the "D" dimension of values is correct:
if (_spaceTime)
{
TEUCHOS_TEST_FOR_EXCEPTION(values.dimension(2) != spaceDim+1, std::invalid_argument, "For space-time normals, values.dimension(2) should be spaceDim + 1.");
}
else
{
TEUCHOS_TEST_FOR_EXCEPTION(values.dimension(2) != spaceDim, std::invalid_argument, "For spatial normals, values.dimension(2) should be spaceDim.");
}
}
const Intrepid::FieldContainer<double> *sideNormals = _spaceTime ? &(basisCache->getSideNormalsSpaceTime()) : &(basisCache->getSideNormals());
int comp = _comp;
if (comp == -2)
{
// want to select the temporal component, t()
comp = spaceDim;
}
for (int cellIndex=0; cellIndex<numCells; cellIndex++)
{
for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
{
if (comp == -1)
{
for (int d=0; d<spaceDim; d++)
{
double nd = (*sideNormals)(cellIndex,ptIndex,d);
values(cellIndex,ptIndex,d) = nd;
}
if (_spaceTime)
{
double nd = (*sideNormals)(cellIndex,ptIndex,spaceDim);
values(cellIndex,ptIndex,spaceDim) = nd;
}
}
else
{
double ni = (*sideNormals)(cellIndex,ptIndex,comp);
values(cellIndex,ptIndex) = ni;
}
}
}
}
示例2: values
void SideParityFunction::values(Intrepid::FieldContainer<double> &values, BasisCachePtr sideBasisCache)
{
this->CHECK_VALUES_RANK(values);
int numCells = values.dimension(0);
int numPoints = values.dimension(1);
int sideIndex = sideBasisCache->getSideIndex();
if (sideIndex == -1)
{
TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "non-sideBasisCache passed into SideParityFunction");
}
if (sideBasisCache->getCellSideParities().size() > 0)
{
// then we'll use this, and won't require that mesh and cellIDs are set
if (sideBasisCache->getCellSideParities().dimension(0) != numCells)
{
TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "sideBasisCache->getCellSideParities() is non-empty, but the cell dimension doesn't match that of the values FieldContainer.");
}
for (int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++)
{
int parity = sideBasisCache->getCellSideParities()(cellOrdinal,sideIndex);
for (int ptOrdinal=0; ptOrdinal<numPoints; ptOrdinal++)
{
values(cellOrdinal,ptOrdinal) = parity;
}
}
}
else
{
vector<GlobalIndexType> cellIDs = sideBasisCache->cellIDs();
if (cellIDs.size() != numCells)
{
TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "cellIDs.size() != numCells");
}
Teuchos::RCP<Mesh> mesh = sideBasisCache->mesh();
if (! mesh.get())
{
TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "mesh unset in BasisCache.");
}
for (int cellIndex=0; cellIndex<numCells; cellIndex++)
{
int parity = mesh->cellSideParitiesForCell(cellIDs[cellIndex])(0,sideIndex);
for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
{
values(cellIndex,ptIndex) = parity;
}
}
}
}
示例3: mapDataIntoGlobalContainer
void SubBasisDofMatrixMapper::mapDataIntoGlobalContainer(const Intrepid::FieldContainer<double> &allLocalData, const vector<int> &basisOrdinalsInLocalData,
const map<GlobalIndexType, unsigned> &globalIndexToOrdinal,
bool fittableDofsOnly, const set<GlobalIndexType> &fittableDofIndices, Intrepid::FieldContainer<double> &globalData)
{
const set<int>* basisOrdinalFilter = &this->basisDofOrdinalFilter();
vector<int> dofIndices(basisOrdinalFilter->begin(),basisOrdinalFilter->end());
FieldContainer<double> subBasisData(basisOrdinalFilter->size());
int dofCount = basisOrdinalFilter->size();
if (allLocalData.rank()==1)
{
for (int i=0; i<dofCount; i++)
{
subBasisData[i] = allLocalData[basisOrdinalsInLocalData[dofIndices[i]]];
}
}
else
{
TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "mapDataIntoGlobalContainer only supports rank 1 basis data");
}
// subBasisData must be rank 2, and must have the same size as FilteredLocalDofOrdinals in its first dimension
// reshape as a rank 2 container (column vector as a matrix):
subBasisData.resize(subBasisData.dimension(0),1);
this->mapSubBasisDataIntoGlobalContainer(subBasisData, globalIndexToOrdinal, fittableDofsOnly, fittableDofIndices, globalData);
}
示例4:
void ConstantScalarFunction<Scalar>::values(Intrepid::FieldContainer<Scalar> &values, BasisCachePtr basisCache)
{
this->CHECK_VALUES_RANK(values);
for (int i=0; i < values.size(); i++)
{
values[i] = _value;
}
}
示例5:
void ConstantVectorFunction<Scalar>::values(Intrepid::FieldContainer<Scalar> &values, BasisCachePtr basisCache)
{
this->CHECK_VALUES_RANK(values);
int spaceDim = values.dimension(2);
if (spaceDim > _value.size())
{
TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "spaceDim is greater than length of vector...");
}
// values are stored in (C,P,D) order, the important thing here being that we can do this:
for (int i=0; i < values.size(); )
{
for (int d=0; d < spaceDim; d++)
{
values[i++] = _value[d];
}
}
}
示例6: conditionNumberLAPACK
double conditionNumberLAPACK(const Epetra_RowMatrix &stiffnessMatrix, bool diagScaling)
{
Intrepid::FieldContainer<double> A;
SerialDenseWrapper::extractFCFromEpetra_RowMatrix(stiffnessMatrix, A);
if (diagScaling)
{
for (int i=0; i<A.dimension(0); i++)
{
double diag = A(i,i);
for (int j=0; j<A.dimension(1); j++)
{
A(i,j) /= diag;
}
}
}
return SerialDenseWrapper::condest(A);
}
示例7: values
void MaxFunction::values(Intrepid::FieldContainer<double> &values, BasisCachePtr basisCache)
{
this->CHECK_VALUES_RANK(values);
Intrepid::FieldContainer<double> values2(values);
_f1->values(values,basisCache);
_f2->values(values2,basisCache);
for(int i = 0; i < values.size(); i++)
{
values[i] = std::max(values[i],values2[i]);
}
}
示例8: getInterpolatoryCoordinates
/** Get the local coordinates for this field. This is independent of element
* locations.
*
* \param[in,out] coords Coordinates associated with this field type.
*/
void IntrepidFieldPattern::getInterpolatoryCoordinates(const Intrepid::FieldContainer<double> & cellVertices,
Intrepid::FieldContainer<double> & coords) const
{
TEUCHOS_ASSERT(cellVertices.rank()==3);
int numCells = cellVertices.dimension(0);
// grab the local coordinates
Intrepid::FieldContainer<double> localCoords;
getInterpolatoryCoordinates(localCoords);
// resize the coordinates field container
coords.resize(numCells,localCoords.dimension(0),getDimension());
if(numCells>0) {
// map to phsyical coordinates
Intrepid::CellTools<double> cellTools;
cellTools.mapToPhysicalFrame(coords,localCoords,cellVertices,intrepidBasis_->getBaseCellTopology());
}
}
示例9: sqrt
void PolarizedFunction<Scalar>::values(Intrepid::FieldContainer<Scalar> &values, BasisCachePtr basisCache)
{
this->CHECK_VALUES_RANK(values);
static const double PI = 3.141592653589793238462;
int numCells = values.dimension(0);
int numPoints = values.dimension(1);
const Intrepid::FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
Intrepid::FieldContainer<double> polarPoints = *points;
for (int cellIndex=0; cellIndex<numCells; cellIndex++)
{
for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
{
double x = (*points)(cellIndex,ptIndex,0);
double y = (*points)(cellIndex,ptIndex,1);
double r = sqrt(x * x + y * y);
double theta = (r != 0) ? acos(x/r) : 0;
// now x = r cos theta, but need to guarantee that y = r sin theta (might differ in sign)
// according to the acos docs, theta will be in [0, pi], so the rule is: (y < 0) ==> theta := 2 pi - theta;
if (y < 0) theta = 2*PI-theta;
polarPoints(cellIndex, ptIndex, 0) = r;
polarPoints(cellIndex, ptIndex, 1) = theta;
// if (r == 0) {
// cout << "r == 0!" << endl;
// }
}
}
BasisCachePtr dummyBasisCache = Teuchos::rcp( new PhysicalPointCache( polarPoints ) );
_f->values(values,dummyBasisCache);
if (_f->isZero())
{
cout << "Warning: in PolarizedFunction, we are being asked for values when _f is zero. This shouldn't happen.\n";
}
// cout << "polarPoints: \n" << polarPoints;
// cout << "PolarizedFunction, values: \n" << values;
}
示例10: fillFieldContainer
void fillFieldContainer(int fieldNum,const std::string & blockId,
const panzer::UniqueGlobalIndexer<short,int> & ugi,
Intrepid::FieldContainer<int> & data)
{
data.resize(1,4);
const std::vector<short> & elements = ugi.getElementBlock(blockId);
const std::vector<int> & fieldOffsets = ugi.getGIDFieldOffsets(blockId,fieldNum);
std::vector<int> gids;
for(std::size_t e=0;e<elements.size();e++) {
ugi.getElementGIDs(elements[e],gids);
for(std::size_t f=0;f<fieldOffsets.size();f++)
data(e,f) = gids[fieldOffsets[f]];
}
}
示例11: if
void SimpleVectorFunction<Scalar>::values(Intrepid::FieldContainer<Scalar> &values, BasisCachePtr basisCache)
{
this->CHECK_VALUES_RANK(values);
int numCells = values.dimension(0);
int numPoints = values.dimension(1);
const Intrepid::FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
int spaceDim = points->dimension(2);
for (int cellIndex=0; cellIndex<numCells; cellIndex++)
{
for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
{
if (spaceDim == 1)
{
double x = (*points)(cellIndex,ptIndex,0);
values(cellIndex,ptIndex,0) = value(x)[0];
}
else if (spaceDim == 2)
{
double x = (*points)(cellIndex,ptIndex,0);
double y = (*points)(cellIndex,ptIndex,1);
values(cellIndex,ptIndex,0) = value(x,y)[0];
values(cellIndex,ptIndex,1) = value(x,y)[1];
}
else if (spaceDim == 3)
{
double x = (*points)(cellIndex,ptIndex,0);
double y = (*points)(cellIndex,ptIndex,1);
double z = (*points)(cellIndex,ptIndex,2);
values(cellIndex,ptIndex,0) = value(x,y,z)[0];
values(cellIndex,ptIndex,1) = value(x,y,z)[1];
values(cellIndex,ptIndex,2) = value(x,y,z)[2];
}
}
}
}
示例12: vals
void ATO::Integrator<T>::getMeasure(
T& measure,
const Intrepid::FieldContainer<T>& topoVals,
const Intrepid::FieldContainer<T>& coordCon,
T zeroVal, C compare)
//******************************************************************************//
{
typedef unsigned int uint;
int nDims = coordCon.dimension(1);
if(nDims == 2){
std::vector< Vector3D > points;
std::vector< Tri > tris;
// if there are topoVals that are exactly equal to or very near zeroVal,
// there will be all sorts of special cases. If necessary, nudge values
// away from zeroVal.
Intrepid::FieldContainer<T> vals(topoVals);
int nvals = vals.dimension(0);
for(int i=0; i<nvals; i++){
if( fabs(vals(i) - zeroVal) < 1e-9 ) vals(i) = zeroVal + 1e-9;
}
// compute surface mesh in physicsl coordinates
getSurfaceTris(points,tris,
vals,coordCon,
zeroVal,compare);
// foreach tri, add measure
int ntris = tris.size();
for(int i=0; i<ntris; i++){
T minc = getTriMeasure(points,tris[i]);
measure += minc;
}
}
}
示例13: newpoint
void ATO::Integrator<T>::getSurfaceTris(
std::vector< Vector3D >& points,
std::vector< Tri >& tris,
const Intrepid::FieldContainer<T>& topoVals,
const Intrepid::FieldContainer<T>& coordCon,
T zeroVal, C compare)
//******************************************************************************//
{
const CellTopologyData& cellData = *(cellTopology->getBaseCellTopologyData());
// find intersections
std::vector<Intersection> intersections;
uint nEdges = cellData.edge_count;
int nDims = coordCon.dimension(1);
for(int edge=0; edge<nEdges; edge++){
uint i = cellData.edge[edge].node[0], j = cellData.edge[edge].node[1];
if((topoVals(i)-zeroVal)*(topoVals(j)-zeroVal) < 0.0){
Vector3D newpoint(Intrepid::ZEROS);
T factor = fabs(topoVals(i)-zeroVal)/(fabs(topoVals(i)-zeroVal)+fabs(topoVals(j)-zeroVal));
for(uint k=0; k<nDims; k++) newpoint(k) = (1.0-factor)*coordCon(i,k) + factor*coordCon(j,k);
std::pair<int,int> newIntx(i,j);
if(topoVals(i) > zeroVal){int tmp=newIntx.first; newIntx.first=newIntx.second; newIntx.second=tmp;}
intersections.push_back(Intersection(newpoint,newIntx));
}
}
std::vector<std::pair<Vector3D,Vector3D> > segment;
// if there are four intersections, then there are two interfaces:
if( intersections.size() == 4 ){
segment.resize(2);
int numNodes = basis->getCardinality();
RealType cntrVal = 0.0;
for(int node=0; node<numNodes; node++)
cntrVal += topoVals(node);
cntrVal /= numNodes;
// if topoVal at centroid is negative, interssected segments share a positive valued node
if( cntrVal < zeroVal ){
int second = intersections[0].connect.second;
if( second == intersections[1].connect.second ){
segment[0].first = intersections[0].point; segment[0].second = intersections[1].point;
segment[1].first = intersections[2].point; segment[1].second = intersections[3].point;
} else
if( second == intersections[2].connect.second ){
segment[0].first = intersections[0].point; segment[0].second = intersections[2].point;
segment[1].first = intersections[1].point; segment[1].second = intersections[3].point;
} else
if( second == intersections[3].connect.second ){
segment[0].first = intersections[0].point; segment[0].second = intersections[3].point;
segment[1].first = intersections[1].point; segment[1].second = intersections[2].point;
}
} else {
int first = intersections[0].connect.first;
if( first == intersections[1].connect.first ){
segment[0].first = intersections[0].point; segment[0].second = intersections[1].point;
segment[1].first = intersections[2].point; segment[1].second = intersections[3].point;
} else
if( first == intersections[2].connect.first ){
segment[0].first = intersections[0].point; segment[0].second = intersections[2].point;
segment[1].first = intersections[1].point; segment[1].second = intersections[3].point;
} else
if( first == intersections[3].connect.first ){
segment[0].first = intersections[0].point; segment[0].second = intersections[3].point;
segment[1].first = intersections[1].point; segment[1].second = intersections[2].point;
}
}
} else
if( intersections.size() == 2){
segment.resize(1);
segment[0].first = intersections[0].point; segment[0].second = intersections[1].point;
}
std::vector< Teuchos::RCP<MiniPoly> > polys;
int npoints = coordCon.dimension(0), ndims = coordCon.dimension(1);
Teuchos::RCP<MiniPoly> poly = Teuchos::rcp(new MiniPoly(npoints));
std::vector<Vector3D>& pnts = poly->points;
std::vector<int>& map = poly->mapToBase;
for(int pt=0; pt<npoints; pt++){
for(int dim=0; dim<ndims; dim++) pnts[pt](dim) = coordCon(pt,dim);
map[pt] = pt;
}
polys.push_back(poly);
int nseg = segment.size();
for(int seg=0; seg<nseg; seg++)
partitionBySegment(polys, segment[seg]);
typename std::vector< Teuchos::RCP<MiniPoly> >::iterator itpoly;
for(itpoly=polys.begin(); itpoly!=polys.end(); itpoly++)
if( included(*itpoly,topoVals,zeroVal,compare) ) trisFromPoly(points, tris, *itpoly);
}
示例14: topoVals
void ATO::Integrator<T>::getCubature(std::vector<std::vector<T> >& refPoints,
std::vector<T>& weights,
const Intrepid::FieldContainer<T>& coordCon,
const Intrepid::FieldContainer<T>& topoVals, T zeroVal)
//******************************************************************************//
{
// check for degeneracy (any topoVals == zeroVal). These will have to be handled specially.
int nTopoVals = topoVals.dimension(0);
T mult = topoVals(0)-zeroVal;
for(int i=1; i<nTopoVals; i++) mult *= (topoVals(i)-zeroVal);
TEUCHOS_TEST_FOR_EXCEPTION(mult == 0.0, std::runtime_error,
std::endl << "Degenerate topology: handling not yet implemented." << std::endl);
// This function computes the weights and refPoints for the base topology,
// i.e., it ignores subcells associated with an extended topology (see Shards
// Doxygen for details on base versus extended topologies). For example,
// the extra subcells associated with a Quad8 or Quad9 element will be
// ignored, so it will be treated as a Quad4 element.
typedef unsigned int uint;
int nDims = coordCon.dimension(1);
if(nDims == 2){
/*
const CellTopologyData& cellData = *(cellTopology->getBaseCellTopologyData());
std::vector< Intrepid::Vector<T> > points;
std::vector< Intrepid::Vector<int> > tris;
getSurfaceTris(points,tris,
cellData,coordCon,topoVals,zeroVal);
std::vector<Intrepid::Vector<T> > Apoints, Bpoints;
// add negative/positive points to Apoints/Bpoints vector
Intrepid::Vector<T> point(nDims);
for(int i=0; i<nTopoVals; i++){
for(int j=0; j<nDims; j++) point(j) = coordCon(i,j);
if(topoVals(i) < zeroVal) Apoints.push_back(point);
else Bpoints.push_back(point);
}
// find/count intersection points
uint nIntersections = 0;
uint nEdges = cellData.edge_count;
for(int edge=0; edge<nEdges; edge++){
uint i = cellData.edge[edge].node[0],
j = cellData.edge[edge].node[1];
if((topoVals(i)-zeroVal)*(topoVals(j)-zeroVal) < 0.0){
Intrepid::Vector<T> newpoint(nDims);
T factor = fabs(topoVals(i))/(fabs(topoVals(i))+fabs(topoVals(j)));
for(uint k=0; k<nDims; k++) newpoint(k) = (1.0-factor)*coordCon(i,k) + factor*coordCon(j,k);
Apoints.push_back(newpoint);
Bpoints.push_back(newpoint);
nIntersections++;
}
}
// if there are four intersections, then there are two interfaces:
if( nIntersections == 4 ){
} else
if( nIntersections == 2 ){
// find centerpoint
Intrepid::Vector<T> Acenter(Apoints[0]), Bcenter(Bpoints[0]);
uint nApoints = Apoints.size(), nBpoints = Bpoints.size();
for(uint i=1; i<nApoints; i++) Acenter += Apoints[i];
Acenter /= nApoints;
for(uint i=1; i<nBpoints; i++) Bcenter += Bpoints[i];
Bcenter /= nBpoints;
// sort by counterclockwise angle about surface normal
T pi = acos(-1.0);
Intrepid::Vector<T> Ax(Apoints[0]-Acenter);
Intrepid::Vector<T> Ay(-Ax(1),Ax(0),0.0);
T Arefnorm = Intrepid::norm(Ax);
std::map<T, uint> Aangles;
Aangles.insert( std::pair<T, uint>(0.0,0) );
for(int i=1; i<nApoints; i++){
Intrepid::Vector<T> comp = Apoints[i] - Acenter;
T compnorm = Intrepid::norm(comp);
T angle = acos(Ax * comp/(compnorm*Arefnorm));
if( Ay * comp < 0.0 ) angle = 2.0*pi - angle;
Aangles.insert( std::pair<T, uint>(angle,i) );
}
// create polygons
nApoints = Aangles.size();
if(nApoints == 5){
Intrepid::Vector<T>& c0 = Acenter;
typename std::map<T,uint>::iterator it=Aangles.begin();
while(it!=Aangles.end()){
Intrepid::Vector<T>& c1 = Apoints[it->second];
it++;
Intrepid::Vector<T>& c2 = Apoints[it->second];
addCubature(refPoints, weights, c0, c1, c2);
}
it=Aangles.begin();
//.........这里部分代码省略.........
示例15: main
//.........这里部分代码省略.........
VarPtr uhat = vf->traceVar("uhat");
VarPtr fhat = vf->fluxVar("fhat");
VarPtr u = vf->fieldVar("u");
VarPtr sigma = vf->fieldVar("sigma", VECTOR_L2);
//////////////////// DEFINE BILINEAR FORM ///////////////////////
BFPtr bf = Teuchos::rcp( new BF(vf) );
// tau terms:
bf->addTerm(sigma, tau);
bf->addTerm(u, tau->div());
bf->addTerm(-uhat, tau->dot_normal());
// v terms:
bf->addTerm( sigma, v->grad() );
bf->addTerm( fhat, v);
//////////////////// DEFINE INNER PRODUCT(S) ///////////////////////
IPPtr ip = bf->graphNorm();
//////////////////// SPECIFY RHS ///////////////////////
RHSPtr rhs = RHS::rhs();
FunctionPtr one = Function::constant(1.0);
rhs->addTerm( one * v );
//////////////////// CREATE BCs ///////////////////////
BCPtr bc = BC::bc();
FunctionPtr zero = Function::zero();
SpatialFilterPtr entireBoundary = SpatialFilter::allSpace();
bc->addDirichlet(uhat, entireBoundary, zero);
//////////////////// SOLVE & REFINE ///////////////////////
// Output solution
Intrepid::FieldContainer<GlobalIndexType> savedCellPartition;
Teuchos::RCP<Epetra_FEVector> savedLHSVector;
{
//////////////////// BUILD MESH ///////////////////////
int H1Order = 4, pToAdd = 2;
Teuchos::RCP<Mesh> mesh = Teuchos::rcp( new Mesh (meshTopology, bf, H1Order, pToAdd) );
Teuchos::RCP<Solution> solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );
solution->solve(false);
RefinementStrategy refinementStrategy( solution, 0.2);
HDF5Exporter exporter(mesh, "Poisson");
// exporter.exportSolution(solution, vf, 0, 2, cellIDToSubdivision(mesh, 4));
exporter.exportSolution(solution, 0, 2);
mesh->saveToHDF5("MeshSave.h5");
solution->saveToHDF5("SolnSave.h5");
solution->save("PoissonProblem");
// int numRefs = 1;
// for (int ref = 1; ref <= numRefs; ref++)
// {
// refinementStrategy.refine(commRank==0);
// solution->solve(false);
// mesh->saveToHDF5("MeshSave.h5");
// solution->saveToHDF5("SolnSave.h5");
// exporter.exportSolution(solution, vf, ref, 2, cellIDToSubdivision(mesh, 4));
// }
mesh->globalDofAssignment()->getPartitions(savedCellPartition);
savedLHSVector = solution->getLHSVector();
}
{
SolutionPtr loadedSolution = Solution::load(bf, "PoissonProblem");
HDF5Exporter exporter(loadedSolution->mesh(), "ProblemLoaded");
// exporter.exportSolution(loadedSolution, vf, 0, 2, cellIDToSubdivision(loadedSolution->mesh(), 4));