Write mesh to output file.
115{
118 "ProcessDisplacement needs bnd parameter with a single id.");
119
120 string toFile =
m_config[
"to"].as<
string>();
121
122 if (toFile == "")
123 {
124 cout << "ProcessDisplacement: you must provide a file" << endl;
125 return;
126 }
127
128 bool useVertexIds =
m_config[
"usevertexids"].as<
bool>();
129
130 vector<string> files;
131 files.push_back(toFile);
136
137
138 int bndCondId =
m_config[
"bnd"].as<
int>();
139
140
141
142 if (bndGraph->GetMeshDimension() == 1)
143 {
144 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"v"));
145 m_f->m_variables.push_back(
"v");
146
148 m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
150 m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
151
152 map<int, int> bndCondIds;
153 for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
154 {
155 bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()] = i;
156 }
157
159
160 for (auto &sIt : tmp)
161 {
162 auto mIt = bndCondIds.find(sIt.first);
163
164 if (mIt == bndCondIds.end())
165 {
166 cout << "Warning: couldn't find element " << sIt.first << endl;
167 continue;
168 }
169
170 int e = mIt->second;
171
173 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
174 bndCondExpU->GetExp(e)->GetGeom());
175
177
178
181 bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(), to);
182
183 const int offset = bndCondExpU->GetPhys_Offset(e);
184 const int nq = toSeg->GetTotPoints();
185
186 Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
187
188 bndCondExpU->GetExp(e)->GetCoords(xC, yC);
189 toSeg->GetCoords(xL, yL);
190
192 tmp = bndCondExpU->UpdatePhys() + offset, 1);
194 tmp = bndCondExpV->UpdatePhys() + offset, 1);
195 }
196
197
198 bndCondExpU->FwdTransBndConstrained(bndCondExpU->GetPhys(),
199 bndCondExpU->UpdateCoeffs());
200 bndCondExpV->FwdTransBndConstrained(bndCondExpV->GetPhys(),
201 bndCondExpV->UpdateCoeffs());
202 }
203 else if (bndGraph->GetMeshDimension() == 2)
204 {
205 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"v"));
206 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"w"));
207 m_f->m_variables.push_back(
"v");
208 m_f->m_variables.push_back(
"w");
209
211 m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
213 m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
215 m_f->m_exp[2]->GetBndCondExpansions()[bndCondId];
216
217 map<int, int> bndCondIds;
218 for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
219 {
220 bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()] = i;
221 }
222
224
225 if (useVertexIds)
226 {
227 for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
228 {
230 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
231 bndCondExpU->GetExp(i)->GetGeom());
232
233 TriFaceIDs t(from->GetVid(0), from->GetVid(1), from->GetVid(2));
234 vertexFaceMap[t] = i;
235 }
236 }
237
239
240 for (auto &sIt : tmp)
241 {
242 int e;
243
244 if (useVertexIds)
245 {
246 TriFaceIDs t(sIt.second->GetVid(0), sIt.second->GetVid(1),
247 sIt.second->GetVid(2));
248
249 auto tIt = vertexFaceMap.find(t);
250 e = tIt == vertexFaceMap.end() ? -1 : tIt->second;
251 }
252 else
253 {
254 auto mIt = bndCondIds.find(sIt.first);
255 e = mIt == bndCondIds.end() ? -1 : mIt->second;
256 }
257
258 if (e == -1)
259 {
260 cout << "Warning: couldn't find element " << sIt.first << endl;
261 continue;
262 }
263
265 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
266 bndCondExpU->GetExp(e)->GetGeom());
267
269
270
273 bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
274 bndCondExpV->GetExp(e)->GetBasis(1)->GetBasisKey(), to);
275
276 const int offset = bndCondExpU->GetPhys_Offset(e);
277 const int nq = toSeg->GetTotPoints();
278
279 Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
280 Array<OneD, NekDouble> zL(nq), zC(nq);
281
282 bndCondExpU->GetExp(e)->GetCoords(xC, yC, zC);
283 toSeg->GetCoords(xL, yL, zL);
284
286 tmp = bndCondExpU->UpdatePhys() + offset, 1);
288 tmp = bndCondExpV->UpdatePhys() + offset, 1);
290 tmp = bndCondExpW->UpdatePhys() + offset, 1);
291 }
292
293
294 bndCondExpU->FwdTransBndConstrained(bndCondExpU->GetPhys(),
295 bndCondExpU->UpdateCoeffs());
296 bndCondExpV->FwdTransBndConstrained(bndCondExpV->GetPhys(),
297 bndCondExpV->UpdateCoeffs());
298 bndCondExpW->FwdTransBndConstrained(bndCondExpW->GetPhys(),
299 bndCondExpW->UpdateCoeffs());
300 }
301}
#define ASSERTL0(condition, msg)
FieldSharedPtr m_f
Field object.
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true, SpatialDomains::MeshGraphSharedPtr partitionedGraph=nullptr)
std::unordered_map< TriFaceIDs, int, TriFaceHash > TriFaceMap
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< SegExp > SegExpSharedPtr
std::shared_ptr< TriExp > TriExpSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::map< int, TriGeomSharedPtr > TriGeomMap
std::map< int, SegGeomSharedPtr > SegGeomMap
std::shared_ptr< SegGeom > SegGeomSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.