44 #include <boost/math/special_functions/fpclassify.hpp>
53 ModuleKey ProcessEquiSpacedOutput::className =
56 ProcessEquiSpacedOutput::create,
57 "Write data as equi-spaced output using simplices to represent the data for connecting points");
63 f->m_setUpEquiSpacedFields =
true;
66 "Only process tetrahedral elements");
87 cout <<
"Interpolating fields to equispaced" << endl;
90 int coordim =
m_f->m_exp[0]->GetCoordim(0);
91 int shapedim =
m_f->m_exp[0]->GetExp(0)->GetShapeDimension();
92 int npts =
m_f->m_exp[0]->GetTotPoints();
95 int nel =
m_f->m_exp[0]->GetExpSize();
98 bool homogeneous1D =
false;
99 if (
m_f->m_fielddef.size())
101 if (
m_f->m_fielddef[0]->m_numHomogeneousDir == 1)
103 ASSERTL0(shapedim==2,
"Homogeneous only implemented for 3DH1D");
106 homogeneous1D =
true;
108 else if(
m_f->m_fielddef[0]->m_numHomogeneousDir == 2)
110 ASSERTL0(
false,
"Homegeneous2D case not supported");
115 if(
m_f->m_session->DefinesSolverInfo(
"HOMOGENEOUS"))
117 std::string HomoStr =
m_f->m_session->GetSolverInfo(
"HOMOGENEOUS");
119 if((HomoStr ==
"HOMOGENEOUS1D") || (HomoStr ==
"Homogeneous1D")
120 || (HomoStr ==
"1D") || (HomoStr ==
"Homo1D"))
122 ASSERTL0(shapedim==2,
"Homogeneous only implemented for 3DH1D");
125 homogeneous1D =
true;
127 if((HomoStr ==
"HOMOGENEOUS2D") || (HomoStr ==
"Homogeneous2D")
128 || (HomoStr ==
"2D") || (HomoStr ==
"Homo2D"))
130 ASSERTL0(
false,
"Homegeneous2D case not supported");
137 int newtotpoints = 0;
145 map<int,StdRegions::Orientation > face0orient;
146 set<int> prismorient;
150 vector<std::string> fieldNames;
152 vector<Array<OneD, int> > ptsConn;
155 for(
int i = 0; i < nel; ++i)
157 e =
m_f->m_exp[0]->GetExp(i);
161 int fid = e->GetGeom()->GetFid(0);
162 if(face0orient.count(fid))
164 prismorient.insert(i);
184 for(
int i = 0; i < nel; ++i)
186 e =
m_f->m_exp[0]->GetExp(i);
189 int fid = e->GetGeom()->GetFid(2);
191 if(face0orient.count(fid))
206 prismorient.insert(i);
213 prismorient.insert(i);
220 for(
int i = 0; i < nel; ++i)
222 e =
m_f->m_exp[0]->GetExp(i);
225 if(
m_f->m_exp[0]->GetExp(i)->DetShapeType() !=
232 switch(e->DetShapeType())
236 int npoints0 = e->GetBasis(0)->GetNumPoints();
244 int np0 = e->GetBasis(0)->GetNumPoints();
245 int np1 = e->GetBasis(1)->GetNumPoints();
246 int np = max(np0,np1);
253 int np0 = e->GetBasis(0)->GetNumPoints();
254 int np1 = e->GetBasis(1)->GetNumPoints();
255 int np = max(np0,np1);
263 int np0 = e->GetBasis(0)->GetNumPoints();
264 int np1 = e->GetBasis(1)->GetNumPoints();
265 int np2 = e->GetBasis(2)->GetNumPoints();
266 int np = max(np0,max(np1,np2));
274 int np0 = e->GetBasis(0)->GetNumPoints();
275 int np1 = e->GetBasis(1)->GetNumPoints();
276 int np2 = e->GetBasis(2)->GetNumPoints();
277 int np = max(np0,max(np1,np2));
285 int np0 = e->GetBasis(0)->GetNumPoints();
286 int np1 = e->GetBasis(1)->GetNumPoints();
287 int np2 = e->GetBasis(2)->GetNumPoints();
288 int np = max(np0,max(np1,np2));
296 int np0 = e->GetBasis(0)->GetNumPoints();
297 int np1 = e->GetBasis(1)->GetNumPoints();
298 int np2 = e->GetBasis(2)->GetNumPoints();
299 int np = max(np0,max(np1,np2));
311 ppe.push_back(newpoints);
312 newtotpoints += newpoints;
317 bool standard =
true;
319 if(prismorient.count(i))
324 e->GetSimplexEquiSpacedConnectivity(conn,standard);
329 if((prevNcoeffs != e->GetNcoeffs()) ||
330 (prevNpoints != e->GetTotPoints()))
332 prevNcoeffs = e->GetNcoeffs();
333 prevNpoints = e->GetTotPoints();
335 e->GetSimplexEquiSpacedConnectivity(conn);
339 for(
int j = 0; j < conn.num_elements(); ++j)
341 newconn[j] = conn[j] + cnt;
344 ptsConn.push_back(newconn);
348 if(
m_f->m_fielddef.size())
350 nfields =
m_f->m_exp.size();
359 for(
int i = 0; i < nfields + coordim; ++i)
365 for(
int i = 0; i < coordim; ++i)
370 for(
int i = coordim; i < 3; ++i)
375 m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
377 int nq1 =
m_f->m_exp[0]->GetTotPoints();
383 m_f->m_exp[0]->GetCoords(x1, y1, z1);
388 for(
int n = 0; n < coordim; ++n)
392 for(
int i = 0; i < nel; ++i)
394 m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
396 tmp = pts[n] + cnt1);
398 cnt +=
m_f->m_exp[0]->GetExp(i)->GetTotPoints();
402 if(
m_f->m_fielddef.size())
404 ASSERTL0(
m_f->m_fielddef[0]->m_fields.size() ==
m_f->m_exp.size(),
405 "More expansion defined than fields");
407 for(
int n = 0; n <
m_f->m_exp.size(); ++n)
412 if(
m_config[
"modalenergy"].m_beenSet)
415 for(
int i = 0; i < nel; ++i)
419 cnt +=
m_f->m_exp[0]->GetExp(i)->GetTotPoints();
425 for(
int i = 0; i < nel; ++i)
427 m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
429 tmp = pts[coordim + n] + cnt1);
431 cnt +=
m_f->m_exp[0]->GetExp(i)->GetTotPoints();
436 fieldNames.push_back(
m_f->m_fielddef[0]->m_fields[n]);
445 else if (shapedim == 3)
449 m_f->m_fieldPts->SetConnectivity(ptsConn);
459 int nel =
m_f->m_exp[0]->GetPlane(0)->GetExpSize();
460 int nPlanes =
m_f->m_exp[0]->GetZIDs().num_elements();
461 int npts =
m_f->m_fieldPts->GetNpoints();
462 int nptsPerPlane = npts/nPlanes;
466 for(
int i = 0; i < nel; ++i)
468 e =
m_f->m_exp[0]->GetPlane(0)->GetExp(i);
470 int np0 = e->GetBasis(0)->GetNumPoints();
471 int np1 = e->GetBasis(1)->GetNumPoints();
472 int np = max(np0,np1);
473 maxN = max(np, maxN);
480 for(
int i = 0; i < nel; ++i)
482 e =
m_f->m_exp[0]->GetPlane(0)->GetExp(i);
483 int np0 = e->GetBasis(0)->GetNumPoints();
484 int np1 = e->GetBasis(1)->GetNumPoints();
485 int np = max(np0,np1);
486 switch(e->DetShapeType())
494 vId[cnt1] = e->GetGeom()->GetVid(0);
495 vId[cnt1+np-1] = e->GetGeom()->GetVid(1);
496 vId[cnt1+newpoints-1] = e->GetGeom()->GetVid(2);
502 for (
int n = 1; n < np-1; n++)
505 edgeOrient = e->GetGeom()->GetEorient(0);
508 vId[cnt1+n] = 4*nel + maxN*e->GetGeom()->GetEid(0) + n;
512 vId[cnt1+np-1-n] = 4*nel +
513 maxN*e->GetGeom()->GetEid(0) + n;
517 edgeOrient = e->GetGeom()->GetEorient(1);
521 vId[cnt1+np-1+edge1] = 4*nel +
522 maxN*e->GetGeom()->GetEid(1) + n;
527 vId[cnt1+newpoints-1-edge1] = 4*nel +
528 maxN*e->GetGeom()->GetEid(1) + n;
532 edgeOrient = e->GetGeom()->GetEorient(2);
536 vId[cnt1+newpoints-1-edge2] = 4*nel +
537 maxN*e->GetGeom()->GetEid(2) + n;
542 vId[cnt1+edge2] = 4*nel +
543 maxN*e->GetGeom()->GetEid(2) + n;
549 for (
int n = 1; n < np-1; n++)
551 for (
int m = 1; m < np-n-1; m++)
553 vId[cnt1+mode] = 4*nel + maxN*4*nel + cnt2;
566 vId[cnt1] = e->GetGeom()->GetVid(0);
567 vId[cnt1+np-1] = e->GetGeom()->GetVid(1);
568 vId[cnt1+np*np-1] = e->GetGeom()->GetVid(2);
569 vId[cnt1+np*(np-1)] = e->GetGeom()->GetVid(3);
573 for (
int n = 1; n < np-1; n++)
576 edgeOrient = e->GetGeom()->GetEorient(0);
579 vId[cnt1+n] = 4*nel + maxN*e->GetGeom()->GetEid(0) + n;
583 vId[cnt1+np-1-n] = 4*nel +
584 maxN*e->GetGeom()->GetEid(0) + n;
588 edgeOrient = e->GetGeom()->GetEorient(1);
591 vId[cnt1+np-1+n*np] = 4*nel +
592 maxN*e->GetGeom()->GetEid(1) + n;
596 vId[cnt1+np*np-1-n*np] = 4*nel +
597 maxN*e->GetGeom()->GetEid(1) + n;
601 edgeOrient = e->GetGeom()->GetEorient(2);
604 vId[cnt1+np*np-1-n] = 4*nel +
605 maxN*e->GetGeom()->GetEid(2) + n;
609 vId[cnt1+np*(np-1)+n] = 4*nel +
610 maxN*e->GetGeom()->GetEid(2) + n;
614 edgeOrient = e->GetGeom()->GetEorient(3);
617 vId[cnt1+np*(np-1)-n*np] = 4*nel +
618 maxN*e->GetGeom()->GetEid(3) + n;
622 vId[cnt1+n*np] = 4*nel +
623 maxN*e->GetGeom()->GetEid(3) + n;
628 for (
int n = 1; n < np-1; n++)
630 for (
int m = 1; m < np-1; m++)
632 vId[cnt1+m*np+n] = 4*nel + maxN*4*nel + cnt2;
647 vector<Array<OneD, int> > oldConn;
648 vector<Array<OneD, int> > newConn;
650 m_f->m_fieldPts->GetConnectivity(oldConn);
652 for(
int i = 0; i < nel; ++i)
654 int nTris = oldConn[i].num_elements()/3;
659 for(
int n = 0; n<nTris; n++)
662 int vId0 = vId[oldConn[i][3*n+0]];
663 int vId1 = vId[oldConn[i][3*n+1]];
664 int vId2 = vId[oldConn[i][3*n+2]];
668 if ( (vId0<vId1) && (vId0<vId2))
682 else if ( (vId1<vId0) && (vId1<vId2))
712 for (
int p = 0; p<nPlanes-1; p++)
714 conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
715 conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[1]];
716 conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[2]];
717 conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
719 conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
720 conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[1]];
721 conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
722 conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
724 conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
725 conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
726 conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
727 conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[0]];
730 newConn.push_back(conn);
732 m_f->m_fieldPts->SetConnectivity(newConn);
742 e =
m_f->m_exp[0]->GetExp(n);
744 switch(e->DetShapeType())
748 int np0 = e->GetBasis(0)->GetNumPoints();
749 int np1 = e->GetBasis(1)->GetNumPoints();
750 int np = max(np0,np1);
764 e->GetBasis(1)->GetBasisKey(),
772 int np0 = e->GetBasis(0)->GetNumPoints();
773 int np1 = e->GetBasis(1)->GetNumPoints();
774 int np = max(np0,np1);
786 e->GetBasis(1)->GetBasisKey(),
793 ASSERTL0(
false,
"Shape needs setting up");
#define ASSERTL0(condition, msg)
int getNumberOfCoefficients(int Na, int Nb, int Nc)
pair< ModuleType, string > ModuleKey
static Array< OneD, NekDouble > NullNekDouble1DArray
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void SetupEquiSpacedField(void)
map< string, ConfigOption > m_config
List of configuration values.
void GenOrthoModes(int n, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &coeffs)
int getNumberOfCoefficients(int Na, int Nb, int Nc)
FieldSharedPtr m_f
Field object.
Principle Orthogonal Functions .
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
int getNumberOfCoefficients(int Na)
int getNumberOfCoefficients(int Na, int Nb)
Principle Orthogonal Functions .
Defines a specification for a set of points.
boost::shared_ptr< Expansion > ExpansionSharedPtr
boost::shared_ptr< Field > FieldSharedPtr
Represents a command-line configuration option.
int getNumberOfCoefficients(int Na, int Nb, int Nc)
virtual ~ProcessEquiSpacedOutput()
int getNumberOfCoefficients(int Na, int Nb)
int getNumberOfCoefficients(int Na, int Nb, int Nc)
void SetHomogeneousConnectivity(void)
Describes the specification for a Basis.
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...