44 namespace MultiRegions
58 m_bndCondExpansions(),
79 const std::string &variable,
80 const bool SetUpJustDG)
82 m_bndCondExpansions(),
133 int ElmtPointGeom = 0;
134 int TracePointGeom = 0;
137 for (
int i = 0; i <
m_exp->size(); ++i)
140 for (
int j = 0; j < exp1d->GetNverts(); ++j)
142 ElmtPointGeom = (exp1d->GetGeom1D())->GetVid(j);
144 for (
int k = 0; k <
m_trace->GetExpSize(); ++k)
146 TracePointGeom =
m_trace->GetExp(k)->GetGeom()->GetVid(0);
148 if (TracePointGeom == ElmtPointGeom)
171 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
176 ASSERTL0(
false,
"Periodic verts need setting up");
185 for (cnt = n = 0; n <
m_exp->size(); ++n)
187 for (
int v = 0; v < (*m_exp)[n]->GetNverts(); ++v, ++cnt)
194 boost::unordered_map<int,pair<int,int> > perVertToExpMap;
195 boost::unordered_map<int,pair<int,int> >
::iterator it2;
196 for (n = 0; n <
m_exp->size(); ++n)
198 for (
int v = 0; v < (*m_exp)[n]->GetNverts(); ++v)
201 (*
m_exp)[n]->GetGeom()->GetVid(v));
205 perVertToExpMap[it->first] = make_pair(n,v);
212 for (n = 0; n <
m_exp->size(); ++n)
214 for (
int v = 0; v < (*m_exp)[n]->GetNverts(); ++v)
216 int vertGeomId = (*m_exp)[n]->GetGeom()->GetVid(v);
224 it2 = perVertToExpMap.find(ent.
id);
226 if (it2 == perVertToExpMap.end())
228 if (
m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
235 ASSERTL1(
false,
"Periodic vert not found!");
239 int offset =
m_trace->GetPhys_Offset((m_traceMap->GetElmtToTrace())
240 [n][v]->GetElmtId());
243 int offset2 =
m_trace->GetPhys_Offset((m_traceMap->GetElmtToTrace())
245 [it2->second.second]->GetElmtId());
261 if (traceEl->GetLeftAdjacentElementVertex () == -1 ||
262 traceEl->GetRightAdjacentElementVertex() == -1)
272 int traceGeomId = traceEl->GetGeom0D()->GetGlobalID();
278 fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
282 int offset =
m_trace->GetPhys_Offset(traceEl->GetElmtId());
284 GetTraceToUniversalMapUnique(offset) >= 0;
288 else if (traceEl->GetLeftAdjacentElementVertex () != -1 &&
289 traceEl->GetRightAdjacentElementVertex() != -1)
293 (traceEl->GetLeftAdjacentElementExp().get()) ==
298 ASSERTL2(
false,
"Unconnected trace element!");
311 const std::string &variable)
318 map<int,int> GeometryToRegionsMap;
320 SpatialDomains::BoundaryRegionCollection::const_iterator it;
328 for(it = bregions.begin(); it != bregions.end(); ++it)
331 for (bregionIt = it->second->begin();
332 bregionIt != it->second->end(); bregionIt++)
336 int id = (*(bregionIt->second))[0]->GetGlobalID();
337 GeometryToRegionsMap[id] = it->first;
342 map<int,SpatialDomains::GeometrySharedPtr> EndOfDomain;
345 for(domIt = domain.begin(); domIt != domain.end(); ++domIt)
348 for(
int i = 0; i < geomvector->size(); ++i)
350 for(
int j = 0; j < 2; ++j)
352 int vid = (*geomvector)[i]->GetVid(j);
353 if(EndOfDomain.count(vid) == 0)
355 EndOfDomain[vid] = (*geomvector)[i]->GetVertex(j);
359 EndOfDomain.erase(vid);
364 ASSERTL1(EndOfDomain.size() == 2,
"Did not find two ends of domain");
368 for(regIt = EndOfDomain.begin(); regIt != EndOfDomain.end(); ++regIt)
370 if(GeometryToRegionsMap.count(regIt->first) != 0)
373 ASSERTL1(iter != GeometryToRegionsMap.end(),
"Failied to find GeometryToRegionMap");
374 int regionId = iter->second;
375 SpatialDomains::BoundaryRegionCollection::const_iterator bregionsIter = bregions.find(regionId);
376 ASSERTL1(bregionsIter != bregions.end(),
"Failed to find boundary region");
378 returnval->AddBoundaryRegions (regionId,breg);
380 SpatialDomains::BoundaryConditionCollection::const_iterator bconditionsIter = bconditions.find(regionId);
381 ASSERTL1(bconditionsIter != bconditions.end(),
"Failed to find boundary collection");
383 returnval->AddBoundaryConditions(regionId,bcond);
391 gvec->push_back(regIt->second);
392 (*breg)[regIt->first] = gvec;
394 returnval->AddBoundaryRegions(bregions.size()+numNewBc,breg);
400 (*bCondition)[variable] = notDefinedCondition;
402 returnval->AddBoundaryConditions(bregions.size()+numNewBc,bCondition);
423 const std::string &variable,
424 bool SetToOneSpaceDimension):
425 ExpList1D(pSession,graph1D,domain, true,variable,SetToOneSpaceDimension),
426 m_bndCondExpansions(),
445 m_bndCondExpansions(In.m_bndCondExpansions),
446 m_bndConditions(In.m_bndConditions),
447 m_globalBndMat(In.m_globalBndMat),
449 m_traceMap(In.m_traceMap),
450 m_boundaryVerts(In.m_boundaryVerts),
451 m_periodicVerts(In.m_periodicVerts),
452 m_periodicFwdCopy(In.m_periodicFwdCopy),
453 m_periodicBwdCopy(In.m_periodicBwdCopy),
454 m_leftAdjacentVerts(In.m_leftAdjacentVerts)
488 const std::string variable)
496 SpatialDomains::BoundaryRegionCollection::const_iterator it;
499 for (it = bregions.begin(); it != bregions.end(); ++it)
503 if (boundaryCondition->GetBoundaryConditionType() !=
507 for (bregionIt = it->second->begin();
508 bregionIt != it->second->end(); bregionIt++)
510 cnt += bregionIt->second->size();
516 = Array<OneD,MultiRegions::ExpListSharedPtr>(cnt);
519 = Array<OneD,SpatialDomains::BoundaryConditionShPtr>(cnt);
537 const std::string variable)
545 = boost::dynamic_pointer_cast<
547 SpatialDomains::BoundaryRegionCollection::const_iterator it;
552 int i, region1ID, region2ID;
556 map<int,int> BregionToVertMap;
560 for (it = bregions.begin(); it != bregions.end(); ++it)
564 if (locBCond->GetBoundaryConditionType()
570 int id = (*(it->second->begin()->second))[0]->GetGlobalID();
572 BregionToVertMap[it->first] = id;
578 int n = vComm->GetSize();
579 int p = vComm->GetRank();
581 Array<OneD, int> nregions(n, 0);
582 nregions[p] = BregionToVertMap.size();
587 Array<OneD, int> regOffset(n, 0);
589 for (i = 1; i < n; ++i)
591 regOffset[i] = regOffset[i-1] + nregions[i-1];
594 Array<OneD, int> bregmap(totRegions, 0);
595 Array<OneD, int> bregid (totRegions, 0);
596 for(i = regOffset[p], iit = BregionToVertMap.begin();
597 iit != BregionToVertMap.end(); ++iit, ++i)
599 bregid [i] = iit->first;
600 bregmap[i] = iit->second;
601 islocal.insert(iit->first);
607 for (
int i = 0; i < totRegions; ++i)
609 BregionToVertMap[bregid[i]] = bregmap[i];
613 for (it = bregions.begin(); it != bregions.end(); ++it)
617 if (locBCond->GetBoundaryConditionType()
624 region1ID = it->first;
625 region2ID = boost::static_pointer_cast<
627 locBCond)->m_connectedBoundaryRegion;
629 ASSERTL0(BregionToVertMap.count(region1ID) != 0,
630 "Cannot determine vertex of region1ID");
632 ASSERTL0(BregionToVertMap.count(region2ID) != 0,
633 "Cannot determine vertex of region2ID");
637 islocal.count(region2ID) != 0);
659 const std::string variable,
660 Array<OneD, MultiRegions::ExpListSharedPtr>
662 Array<
OneD, SpatialDomains
672 SpatialDomains::BoundaryRegionCollection::const_iterator it;
680 for (it = bregions.begin(); it != bregions.end(); ++it)
683 bconditions, it->first, variable);
685 if (locBCond->GetBoundaryConditionType() ==
690 for (bregionIt = it->second->begin();
691 bregionIt != it->second->end(); bregionIt++)
693 for (k = 0; k < bregionIt->second->size(); k++)
695 if ((vert = boost::dynamic_pointer_cast
697 (*bregionIt->second)[k])))
702 bndCondExpansions[cnt] = locPointExp;
703 bndConditions[cnt++] = locBCond;
708 "dynamic cast to a vertex failed");
716 for (it = bregions.begin(); it != bregions.end(); ++it)
720 switch(locBCond->GetBoundaryConditionType())
727 for (bregionIt = it->second->begin();
728 bregionIt != it->second->end(); bregionIt++)
730 for (k = 0; k < bregionIt->second->size(); k++)
732 if((vert = boost::dynamic_pointer_cast
734 (*bregionIt->second)[k])))
739 bndCondExpansions[cnt] = locPointExp;
740 bndConditions[cnt++] = locBCond;
745 "dynamic cast to a vertex failed");
755 ASSERTL0(
false,
"This type of BC not implemented yet");
768 "Routine currently only tested for HybridDGHelmholtz");
771 "Full matrix global systems are not supported for HDG "
776 "The local to global map is not set up for the requested "
785 (*m_globalBndMat)[mkey] = glo_matrix;
789 glo_matrix = matrixIter->second;
800 Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
808 for(
int v = 0; v < 2; ++v)
814 if(vertExp->GetLeftAdjacentElementExp()->GetGeom()->GetGlobalID() != (*m_exp)[i]->GetGeom()->GetGlobalID())
836 Array<OneD, NekDouble> &Bwd)
869 const Array<OneD, const NekDouble> &field,
870 Array<OneD, NekDouble> &Fwd,
871 Array<OneD, NekDouble> &Bwd)
885 int nLocalSolutionPts;
890 Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
903 for (cnt = n = 0; n < nElements; ++n)
911 nLocalSolutionPts = (*m_exp)[n]->GetNumPoints(0);
913 Basis = (*m_exp)[n]->GetBasis(0);
915 for(v = 0; v < 2; ++v, ++cnt)
917 int offset =
m_trace->GetPhys_Offset(elmtToTrace[n][v]->GetElmtId());
921 (*m_exp)[n]->GetVertexPhysVals(v, field + phys_offset,
926 (*m_exp)[n]->GetVertexPhysVals(v, field + phys_offset,
950 "Method not set up for non-zero Neumann "
951 "boundary condition");
966 "Method not set up for this boundary condition.");
984 Array<OneD, NekDouble> &outarray)
1004 const Array<OneD, const NekDouble> &inarray,
1005 Array<OneD, NekDouble> &outarray)
1009 int n,p,offset,phys_offset;
1012 "input array is of insufficient length");
1014 for (n = 0; n < nexp; ++n)
1018 for (p = 0; p < (*m_exp)[n]->GetNverts(); ++p)
1020 offset =
m_trace->GetPhys_Offset(
1021 (
m_traceMap->GetElmtToTrace())[n][p]->GetElmtId());
1022 (*m_exp)[n]->GetVertexPhysVals(p, inarray + phys_offset,
1029 const Array<OneD, const NekDouble> &Fn,
1030 Array<OneD, NekDouble> &outarray)
1032 int n,offset, t_offset;
1034 Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
1045 Basis = (*m_exp)[n]->GetBasis(0);
1048 int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1055 t_offset =
GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
1056 if(negatedFluxNormal[2*n])
1058 outarray[offset] -= Fn[t_offset];
1062 outarray[offset] += Fn[t_offset];
1065 t_offset =
GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
1067 if(negatedFluxNormal[2*n+1])
1069 outarray[offset+(*m_exp)[n]->GetVertexMap(1)] -= Fn[t_offset];
1073 outarray[offset+(*m_exp)[n]->GetVertexMap(1)] += Fn[t_offset];
1089 Array<OneD, NekDouble> coords(3, 0.0);
1093 for(p = 0; p < 2; ++p)
1096 for (
int i=0; i<((*m_exp)[n]->
1097 GetVertexNormal(p)).num_elements(); i++)
1099 vertnorm += ((*m_exp)[n]->GetVertexNormal(p))[i][0];
1100 coords[0] = vertnorm ;
1103 t_offset =
GetTrace()->GetPhys_Offset(n+p);
1105 if (vertnorm >= 0.0)
1107 m_Ixm = BASE->GetI(coords);
1110 for (j = 0; j < e_ncoeffs; j++)
1112 outarray[offset + j] +=
1113 (m_Ixm->GetPtr())[j] * Fn[t_offset];
1119 m_Ixm = BASE->GetI(coords);
1121 for (j = 0; j < e_ncoeffs; j++)
1123 outarray[offset + j] -=
1124 (m_Ixm->GetPtr())[j] * Fn[t_offset];
1131 static int sav_ncoeffs = 0;
1132 if(!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
1142 Array<OneD, NekDouble> coords(1, 0.0);
1145 m_Ixm = BASE->GetI(coords);
1148 m_Ixp = BASE->GetI(coords);
1150 sav_ncoeffs = e_ncoeffs;
1153 t_offset =
GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
1154 if(negatedFluxNormal[2*n])
1156 for (j = 0; j < e_ncoeffs; j++)
1158 outarray[offset + j] -=
1159 (m_Ixm->GetPtr())[j] * Fn[t_offset];
1164 for (j = 0; j < e_ncoeffs; j++)
1166 outarray[offset + j] +=
1167 (m_Ixm->GetPtr())[j] * Fn[t_offset];
1171 t_offset =
GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
1172 if (negatedFluxNormal[2*n+1])
1174 for (j = 0; j < e_ncoeffs; j++)
1176 outarray[offset + j] -=
1177 (m_Ixp->GetPtr())[j] * Fn[t_offset];
1182 for (j = 0; j < e_ncoeffs; j++)
1184 outarray[offset + j] +=
1185 (m_Ixp->GetPtr())[j] * Fn[t_offset];
1195 const Array<OneD, const NekDouble> &inarray,
1196 Array<OneD, NekDouble> &outarray,
1200 const Array<OneD, const NekDouble> &dirForcing)
1206 Array<OneD,NekDouble> e_f, e_l;
1217 int GloBndDofs =
m_traceMap->GetNumGlobalBndCoeffs();
1218 int NumDirBCs =
m_traceMap->GetNumLocalDirBndCoeffs();
1231 Array<OneD,NekDouble> BndSol = Array<OneD,NekDouble>
1235 Array<OneD,NekDouble> BndRhs(GloBndDofs,0.0);
1239 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
1240 Array<OneD, NekDouble> loc_lambda(LocBndCoeffs);
1247 for (cnt = n = 0; n < nexp; ++n)
1249 nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
1251 e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1253 e_l = loc_lambda + cnt;
1258 Floc =
Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
1273 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i);
1278 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i);
1286 if (GloBndDofs - NumDirBCs > 0)
1307 m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
1310 out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
1324 const std::string varName,
1330 Array<OneD, NekDouble> x0(1);
1331 Array<OneD, NekDouble> x1(1);
1332 Array<OneD, NekDouble> x2(1);
1356 (boost::static_pointer_cast<SpatialDomains
1358 ->m_dirichletCondition).Evaluate(x0[0],x1[0],x2[0],time));
1365 (boost::static_pointer_cast<SpatialDomains
1367 ->m_neumannCondition).Evaluate(x0[0],x1[0],x2[0],time));
1373 (boost::static_pointer_cast<SpatialDomains
1375 ->m_robinFunction).Evaluate(x0[0],x1[0],x2[0],time));
1378 (boost::static_pointer_cast<SpatialDomains
1380 ->m_robinPrimitiveCoeff).Evaluate(x0[0],x1[0],x2[0],time));
1389 ASSERTL0(
false,
"This type of BC not implemented yet");
1398 Array<OneD, int> &ElmtID, Array<OneD,int> &VertID)
1400 map<int, int> VertGID;
1406 if (ElmtID.num_elements() != nbcs)
1408 ElmtID = Array<OneD, int>(nbcs,-1);
1412 fill(ElmtID.get(), ElmtID.get()+nbcs, -1);
1415 if (VertID.num_elements() != nbcs)
1417 VertID = Array<OneD, int>(nbcs);
1424 VertGID[Vid] = cnt++;
1432 for(i = 0; i < exp->GetNverts(); ++i)
1434 id = exp->GetGeom()->GetVid(i);
1436 if (VertGID.count(
id) > 0)
1438 bid = VertGID.find(
id)->second;
1439 ASSERTL1(ElmtID[bid] == -1,
"Edge already set");
1447 ASSERTL1(cnt == nbcs,
"Failed to visit all boundary condtiions");
1467 map<int, RobinBCInfoSharedPtr> returnval;
1468 Array<OneD, int> ElmtID,VertID;
1479 Array<OneD, NekDouble> Array_tmp;
1486 VertID[i],Array_tmp = locExpList->GetPhys());
1491 if(returnval.count(elmtid) != 0)
1493 rInfo->next = returnval.find(elmtid)->second;
1495 returnval[elmtid] = rInfo;