50 namespace MultiRegions
77 ExpList1D::ExpList1D():
88 ExpList(In,DeclareCoeffPhysArrays)
97 const std::vector<unsigned int> &eIDs,
98 const bool DeclareCoeffPhysArrays):
99 ExpList(In, eIDs, DeclareCoeffPhysArrays)
143 = graph1D->GetExpansions();
147 SpatialDomains::ExpansionMap::const_iterator expIt;
148 for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
150 if ((SegmentGeom = boost::dynamic_pointer_cast<
152 expIt->second->m_geomShPtr)))
156 seg->SetElmtId(
id++);
157 (*m_exp).push_back(seg);
161 ASSERTL0(
false,
"dynamic cast to a SegGeom failed");
203 const bool DeclareCoeffPhysArrays):
214 = graph1D->GetExpansions();
217 SpatialDomains::ExpansionMap::const_iterator expIt;
218 for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
223 if ((SegmentGeom = boost::dynamic_pointer_cast<
225 expIt->second->m_geomShPtr)))
231 seg->SetElmtId(
id++);
234 (*m_exp).push_back(seg);
238 ASSERTL0(
false,
"dynamic cast to a SegGeom failed");
250 if(DeclareCoeffPhysArrays)
287 const bool DeclareCoeffPhysArrays,
288 const std::string var,
289 bool SetToOneSpaceDimension):
296 SpatialDomains::CompositeMap::const_iterator compIt;
300 = graph1D->GetExpansions(var);
303 for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
305 comp = compIt->second;
308 for(j = 0; j < compIt->second->size(); ++j)
310 SpatialDomains::ExpansionMap::const_iterator expIt;
312 if((SegmentGeom = boost::dynamic_pointer_cast<
314 (*compIt->second)[j])))
317 if((expIt = expansions.find(SegmentGeom->GetGlobalID())) != expansions.end())
321 if(SetToOneSpaceDimension)
324 SegmentGeom->GenerateOneSpaceDimGeom();
337 ASSERTL0(
false,
"Failed to find basis key");
343 ASSERTL0(
false,
"Failed to dynamic cast geometry to SegGeom");
348 seg->SetElmtId(
id++);
351 (*m_exp).push_back(seg);
363 if(DeclareCoeffPhysArrays)
391 const bool DeclareCoeffPhysArrays,
392 const std::string variable):
401 SpatialDomains::CompositeMap::const_iterator compIt;
406 for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
408 comp = compIt->second;
410 for(j = 0; j < compIt->second->size(); ++j)
412 if((SegmentGeom = boost::dynamic_pointer_cast<
414 (*compIt->second)[j])))
424 seg->SetElmtId(
id++);
425 (*m_exp).push_back(seg);
429 ASSERTL0(
false,
"dynamic cast to a SegGeom failed");
444 if(DeclareCoeffPhysArrays)
476 const bool DeclareCoeffPhysArrays,
477 const std::string variable):
480 int i, j, id, elmtid = 0;
492 map<int,int> EdgeDone;
493 map<int,int> NormalSet;
499 for(i = 0; i < bndCond.num_elements(); ++i)
501 if(bndCond[i]->GetBoundaryConditionType()
504 for(j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
507 ->GetExp(j)->GetBasis(0)->GetBasisKey();
508 exp1D = bndConstraint[i]->GetExp(j)->
509 as<LocalRegions::Expansion1D>();
510 segGeom = exp1D->GetGeom1D();
514 edgesDone.insert(segGeom->GetEid());
516 seg->SetElmtId(elmtid++);
517 (*m_exp).push_back(seg);
527 for(i = 0; i < locexp.size(); ++i)
531 for(j = 0; j < locexp[i]->GetNedges(); ++j)
533 segGeom = exp2D->
GetGeom2D()->GetEdge(j);
534 id = segGeom->GetEid();
536 if (edgesDone.count(
id) != 0)
541 it = edgeOrders.find(
id);
543 if (it == edgeOrders.end())
545 edgeOrders.insert(std::make_pair(
id, std::make_pair(
546 segGeom, locexp[i]->DetEdgeBasisKey(j))));
551 = locexp[i]->DetEdgeBasisKey(j);
560 if (np2 >= np1 && nm2 >= nm1)
564 else if (np2 < np1 && nm2 < nm1)
566 it->second.second = edge;
571 "inappropriate number of points/modes (max "
572 "num of points is not set with max order)");
579 pSession->GetComm()->GetRowComm();
580 int nproc = vComm->GetSize();
581 int edgepr = vComm->GetRank();
588 for(i = 0; i < locexp.size(); ++i)
590 eCnt += locexp[i]->GetNedges();
596 edgesCnt[edgepr] = eCnt;
602 for (i = 1; i < nproc; ++i)
604 eTotOffsets[i] = eTotOffsets[i-1] + edgesCnt[i-1];
612 int cntr = eTotOffsets[edgepr];
614 for(i = 0; i < locexp.size(); ++i)
620 for(j = 0; j < nedges; ++j, ++cntr)
623 locexp[i]->DetEdgeBasisKey(j);
624 EdgesTotID [cntr] = exp2D->GetGeom2D()->GetEid(j);
634 for (i = 0; i < totEdgeCnt; ++i)
636 it = edgeOrders.find(EdgesTotID[i]);
638 if (it == edgeOrders.end())
656 if (np2 >= np1 && nm2 >= nm1)
660 else if (np2 < np1 && nm2 < nm1)
662 it->second.second = edge;
667 "inappropriate number of points/modes (max "
668 "num of points is not set with max order)");
673 for (it = edgeOrders.begin(); it != edgeOrders.end(); ++it)
677 seg->SetElmtId(elmtid++);
678 (*m_exp).push_back(seg);
690 if(DeclareCoeffPhysArrays)
720 for(i = 0; i <
m_exp->size(); ++i)
724 m_offset_elmt_id[i] = i;
725 m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
726 m_npoints += (*m_exp)[i]->GetTotPoints();
760 int quad_npoints = elmExp->GetTotPoints();
762 elmExp->GetPointsType(0));
769 int kernel_width = kernel->GetKernelWidth();
776 for(i = 0; i < inarray.num_elements(); i++)
779 kernel->MoveKernelCenter(inarray[i],local_kernel_breaks);
783 kernel->FindMeshUnderKernel(local_kernel_breaks,h,mesh_breaks);
786 int total_nbreaks = local_kernel_breaks.num_elements() +
787 mesh_breaks.num_elements();
790 kernel->Sort(local_kernel_breaks,mesh_breaks,total_breaks);
795 for(j = 0; j < total_breaks.num_elements()-1; j++)
801 for(r = 0; r < quad_points.num_elements(); r++)
803 mapped_quad_points[r]
804 = (quad_points[r] + 1.0) * 0.5 * (b - a) + a;
813 elmExp->GetBasisNumModes(0),u_value);
817 kernel->EvaluateKernel(mapped_quad_points,h,k_value);
820 for(r = 0; r < quad_npoints; r++)
822 integral_value += (b - a) * 0.5 * k_value[r]
823 * u_value[r] * quad_weights[r];
826 outarray[i] = integral_value/h;
857 for(i = 0; i < outarray.num_elements(); i++)
863 int x_size = inarray2.num_elements();
868 for(i = 0; i < x_size; i++ )
870 x_elm[i] = (int)floor(inarray2[i]/h);
874 for(i = 0; i < x_size; i++)
880 while(x_elm[i] >= num_elm)
882 x_elm[i] -= num_elm ;
887 for(i = 0; i < x_size; i++)
889 x_values_cp[i] = (inarray2[i]/h - floor(inarray2[i]/h))*2 - 1.0;
896 for(i = 0; i < nmodes; i++)
899 jacobi_poly.get()+i*x_size,NULL,i,0.0,0.0);
903 for(r = 0; r < nmodes; r++)
905 for(j = 0; j < x_size; j++)
907 int index = ((x_elm[j])*nmodes)+r;
908 outarray[j] += inarray1[index]*jacobi_poly[r][j];
952 for (i = 0; i <
m_exp->size(); ++i)
954 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
956 (*m_exp)[i]->ComputeVertexNormal(j);
977 int i,j,k,e_npoints,offset;
984 ASSERTL1(Vec.num_elements() >= coordim,
985 "Input vector does not have sufficient dimensions to "
989 for(i = 0; i <
m_exp->size(); ++i)
992 e_npoints = (*m_exp)[i]->GetNumPoints(0);
993 normals = (*m_exp)[i]->GetPhysNormals();
999 for(j = 0; j < e_npoints; ++j)
1003 for(k = 0; k < coordim; ++k)
1005 Vn += Vec[k][offset+j]*normals[k*e_npoints + j];
1011 Upwind[offset + j] = Fwd[offset + j];
1015 Upwind[offset + j] = Bwd[offset + j];
1040 int i,j,e_npoints,offset;
1044 for(i = 0; i <
m_exp->size(); ++i)
1047 e_npoints = (*m_exp)[i]->GetNumPoints(0);
1051 for(j = 0; j < e_npoints; ++j)
1054 if(Vn[offset + j] > 0.0)
1056 Upwind[offset + j] = Fwd[offset + j];
1060 Upwind[offset + j] = Bwd[offset + j];
1077 int i,j,k,e_npoints,offset;
1085 ASSERTL1(normals.num_elements() >= coordim,
1086 "Output vector does not have sufficient dimensions to "
1089 for (i = 0; i <
m_exp->size(); ++i)
1096 int edgeNumber = loc_exp->GetLeftAdjacentElementEdge();
1099 e_npoints = (*m_exp)[i]->GetNumPoints(0);
1101 locnormals = loc_elmt->GetEdgeNormal(edgeNumber);
1102 int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
1103 int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
1105 if (e_nmodes != loc_nmodes)
1107 if (loc_exp->GetRightAdjacentElementEdge() >= 0)
1110 loc_exp->GetRightAdjacentElementExp();
1112 int EdgeNumber = loc_exp->GetRightAdjacentElementEdge();
1115 locnormals = loc_elmt->GetEdgeNormal(EdgeNumber);
1120 for (j = 0; j < e_npoints; ++j)
1124 for (k = 0; k < coordim; ++k)
1126 normals[k][offset + j] = -locnormals[k][j];
1135 for (
int p = 0; p < coordim; ++p)
1139 loc_exp->GetBasis(0)->GetPointsKey();
1141 loc_elmt->GetBasis(0)->GetPointsKey();
1151 for (j = 0; j < e_npoints; ++j)
1155 for (k = 0; k < coordim; ++k)
1157 normals[k][offset + j] = normal[k][j];
1168 for (j = 0; j < e_npoints; ++j)
1172 for (k = 0; k < coordim; ++k)
1174 normals[k][offset + j] = locnormals[k][j];
1206 int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
1208 int ntotminus = (nquad0-1);
1214 (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);
1216 outfile <<
" <Piece NumberOfPoints=\""
1217 << ntot <<
"\" NumberOfCells=\""
1218 << ntotminus <<
"\">" << endl;
1219 outfile <<
" <Points>" << endl;
1220 outfile <<
" <DataArray type=\"Float32\" "
1221 <<
"NumberOfComponents=\"3\" format=\"ascii\">" << endl;
1223 for (i = 0; i < ntot; ++i)
1225 for (j = 0; j < 3; ++j)
1227 outfile << setprecision(8) << scientific
1228 << (float)coords[j][i] <<
" ";
1233 outfile <<
" </DataArray>" << endl;
1234 outfile <<
" </Points>" << endl;
1235 outfile <<
" <Cells>" << endl;
1236 outfile <<
" <DataArray type=\"Int32\" "
1237 <<
"Name=\"connectivity\" format=\"ascii\">" << endl;
1238 for (i = 0; i < nquad0-1; ++i)
1240 outfile << i <<
" " << i+1 << endl;
1243 outfile <<
" </DataArray>" << endl;
1244 outfile <<
" <DataArray type=\"Int32\" "
1245 <<
"Name=\"offsets\" format=\"ascii\">" << endl;
1246 for (i = 0; i < ntotminus; ++i)
1248 outfile << i*2+2 <<
" ";
1251 outfile <<
" </DataArray>" << endl;
1252 outfile <<
" <DataArray type=\"UInt8\" "
1253 <<
"Name=\"types\" format=\"ascii\">" << endl;
1254 for (i = 0; i < ntotminus; ++i)
1259 outfile <<
" </DataArray>" << endl;
1260 outfile <<
" </Cells>" << endl;
1261 outfile <<
" <PointData>" << endl;
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
ExpList1D()
The default constructor.
#define ASSERTL0(condition, msg)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
NekOptimize::GlobalOptParamSharedPtr m_globalOptParam
BasisType GetBasisType() const
Return type of expansion basis.
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
PointsType GetPointsType() const
Return type of quadrature.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
int GetExpSize(void)
This function returns the number of elements in the expansion.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual void v_SetUpPhysNormals()
Set up the normals on each expansion.
std::vector< ExpansionSharedPtr > ExpansionVector
Expansion2DSharedPtr GetLeftAdjacentElementExp() const
void v_Upwind(const Array< OneD, const Array< OneD, NekDouble > > &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Upwind the Fwd and Bwd states based on the velocity field given by Vec.
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
boost::shared_ptr< SegExp > SegExpSharedPtr
Base class for all multi-elemental spectral/hp expansions.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
void PeriodicEval(Array< OneD, NekDouble > &inarray1, Array< OneD, NekDouble > &inarray2, NekDouble h, int nmodes, Array< OneD, NekDouble > &outarray)
Evaluates the global spectral/hp expansion at some arbitray set of points.
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Array< OneD, int > m_offset_elmt_id
Array containing the element id m_offset_elmt_id[n] that the n^th consecutive block of data in m_coef...
int GetNumPoints() const
Return points order at which basis is defined.
boost::shared_ptr< Kernel > KernelSharedPtr
PointsManagerT & PointsManager(void)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Defines a specification for a set of points.
boost::shared_ptr< GeometryVector > Composite
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
std::map< int, Composite > CompositeMap
void SetCoeffPhysOffsets(void)
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void v_GetNormals(Array< OneD, Array< OneD, NekDouble > > &normals)
Populate normals with the normals of all expansions.
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
This class is the abstraction of a one-dimensional multi-elemental expansions which is merely a colle...
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
void PostProcess(LibUtilities::KernelSharedPtr kernel, Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray, NekDouble h, int elmId=0)
Performs the post-processing on a specified element.
virtual void v_ReadGlobalOptimizationParameters()
void ReadGlobalOptimizationParameters()
void CreateCollections(Collections::ImplementationType ImpType=Collections::eNoImpType)
Construct collections of elements containing a single element type and polynomial order from the list...
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
int GetNumModes() const
Returns the order of the basis.
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
const StdRegions::StdExpansionVector &locexp);
void SetExpType(ExpansionType Type)
Returns the type of the expansion.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
int GetNedges() const
This function returns the number of edges of the expansion domain.
virtual ~ExpList1D()
Destructor.
Describes the specification for a Basis.
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
std::map< int, ExpansionShPtr > ExpansionMap
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr