Write mesh to output file.
116{
119 "ProcessDisplacement needs bnd parameter with a single id.");
120
121 string toFile =
m_config[
"to"].as<
string>();
122
123 if (toFile == "")
124 {
125 cout << "ProcessDisplacement: you must provide a file" << endl;
126 return;
127 }
128
129 bool useVertexIds =
m_config[
"usevertexids"].as<
bool>();
130
131 vector<string> files;
132 files.push_back(toFile);
137
138
139 int bndCondId =
m_config[
"bnd"].as<
int>();
140
141
142
143 if (bndGraph->GetMeshDimension() == 1)
144 {
145 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"v"));
146 m_f->m_variables.push_back(
"v");
147
149 m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
151 m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
152
153 map<int, int> bndCondIds;
154 for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
155 {
156 bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()] = i;
157 }
158
160
161 for (auto &sIt : tmp)
162 {
163 auto mIt = bndCondIds.find(sIt.first);
164
165 if (mIt == bndCondIds.end())
166 {
167 cout << "Warning: couldn't find element " << sIt.first << endl;
168 continue;
169 }
170
171 int e = mIt->second;
172
174 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
175 bndCondExpU->GetExp(e)->GetGeom());
176
178
179
182 bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(), to);
183
184 const int offset = bndCondExpU->GetPhys_Offset(e);
185 const int nq = toSeg->GetTotPoints();
186
187 Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
188
189 bndCondExpU->GetExp(e)->GetCoords(xC, yC);
190 toSeg->GetCoords(xL, yL);
191
193 tmp = bndCondExpU->UpdatePhys() + offset, 1);
195 tmp = bndCondExpV->UpdatePhys() + offset, 1);
196 }
197
198
199 bndCondExpU->FwdTransBndConstrained(bndCondExpU->GetPhys(),
200 bndCondExpU->UpdateCoeffs());
201 bndCondExpV->FwdTransBndConstrained(bndCondExpV->GetPhys(),
202 bndCondExpV->UpdateCoeffs());
203 }
204 else if (bndGraph->GetMeshDimension() == 2)
205 {
206 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"v"));
207 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"w"));
208 m_f->m_variables.push_back(
"v");
209 m_f->m_variables.push_back(
"w");
210
212 m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
214 m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
216 m_f->m_exp[2]->GetBndCondExpansions()[bndCondId];
217
218 map<int, int> bndCondIds;
219 for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
220 {
221 bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()] = i;
222 }
223
225
226 if (useVertexIds)
227 {
228 for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
229 {
231 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
232 bndCondExpU->GetExp(i)->GetGeom());
233
234 TriFaceIDs t(from->GetVid(0), from->GetVid(1), from->GetVid(2));
235 vertexFaceMap[t] = i;
236 }
237 }
238
240
241 for (auto &sIt : tmp)
242 {
243 int e;
244
245 if (useVertexIds)
246 {
247 TriFaceIDs t(sIt.second->GetVid(0), sIt.second->GetVid(1),
248 sIt.second->GetVid(2));
249
250 auto tIt = vertexFaceMap.find(t);
251 e = tIt == vertexFaceMap.end() ? -1 : tIt->second;
252 }
253 else
254 {
255 auto mIt = bndCondIds.find(sIt.first);
256 e = mIt == bndCondIds.end() ? -1 : mIt->second;
257 }
258
259 if (e == -1)
260 {
261 cout << "Warning: couldn't find element " << sIt.first << endl;
262 continue;
263 }
264
266 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
267 bndCondExpU->GetExp(e)->GetGeom());
268
270
271
274 bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
275 bndCondExpV->GetExp(e)->GetBasis(1)->GetBasisKey(), to);
276
277 const int offset = bndCondExpU->GetPhys_Offset(e);
278 const int nq = toSeg->GetTotPoints();
279
280 Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
281 Array<OneD, NekDouble> zL(nq), zC(nq);
282
283 bndCondExpU->GetExp(e)->GetCoords(xC, yC, zC);
284 toSeg->GetCoords(xL, yL, zL);
285
287 tmp = bndCondExpU->UpdatePhys() + offset, 1);
289 tmp = bndCondExpV->UpdatePhys() + offset, 1);
291 tmp = bndCondExpW->UpdatePhys() + offset, 1);
292 }
293
294
295 bndCondExpU->FwdTransBndConstrained(bndCondExpU->GetPhys(),
296 bndCondExpU->UpdateCoeffs());
297 bndCondExpV->FwdTransBndConstrained(bndCondExpV->GetPhys(),
298 bndCondExpV->UpdateCoeffs());
299 bndCondExpW->FwdTransBndConstrained(bndCondExpW->GetPhys(),
300 bndCondExpW->UpdateCoeffs());
301 }
302}
#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.