Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessDisplacement.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessDisplacement.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Deforms a mesh given input field defining displacement.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 #include <iostream>
38 using namespace std;
39 
40 #include "ProcessDisplacement.h"
41 
42 #include <StdRegions/StdSegExp.h>
43 #include <StdRegions/StdQuadExp.h>
44 #include <StdRegions/StdTriExp.h>
45 #include <LocalRegions/SegExp.h>
46 #include <LocalRegions/TriExp.h>
49 
50 #include <boost/unordered_set.hpp>
51 #include <boost/tuple/tuple.hpp>
52 #include <boost/tuple/tuple_comparison.hpp>
53 
54 namespace Nektar
55 {
56 namespace Utilities
57 {
58  struct TriFaceIDs
59  {
60  TriFaceIDs(int a, int b, int c) : a(a), b(b), c(c) {}
61  int a;
62  int b;
63  int c;
64  };
65 
66  struct TriFaceHash : std::unary_function<TriFaceIDs, std::size_t>
67  {
68  std::size_t operator()(TriFaceIDs const& p) const
69  {
70  std::vector<int> ids(3);
71 
72  ids[0] = p.a;
73  ids[1] = p.b;
74  ids[2] = p.c;
75 
76  std::sort(ids.begin(), ids.end());
77  return boost::hash_range(ids.begin(), ids.end());
78  }
79  };
80 
81  bool operator==(TriFaceIDs const &p1, TriFaceIDs const &p2)
82  {
83  std::vector<int> ids1(3), ids2(3);
84 
85  ids1[0] = p1.a;
86  ids1[1] = p1.b;
87  ids1[2] = p1.c;
88  ids2[0] = p2.a;
89  ids2[1] = p2.b;
90  ids2[2] = p2.c;
91 
92  std::sort(ids1.begin(), ids1.end());
93  std::sort(ids2.begin(), ids2.end());
94 
95  return ids1[0] == ids2[0] &&
96  ids1[1] == ids2[1] &&
97  ids1[2] == ids2[2];
98  }
99 
100  typedef boost::unordered_map<TriFaceIDs, int, TriFaceHash> TriFaceMap;
101 
102  ModuleKey ProcessDisplacement::className =
104  ModuleKey(eProcessModule, "displacement"), ProcessDisplacement::create,
105  "Deform a mesh given an input field defining displacement");
106 
107  ProcessDisplacement::ProcessDisplacement(FieldSharedPtr f) :
108  ProcessModule(f)
109  {
110  m_config["to"] = ConfigOption(
111  false, "", "Name of file containing high order boundary");
112  m_config["id"] = ConfigOption(
113  false, "", "Boundary ID to calculate displacement for");
114  m_config["usevertexids"] = ConfigOption(
115  false, "0", "Use vertex IDs instead of face IDs for matching");
116  f->m_declareExpansionAsContField = true;
117  f->m_writeBndFld = true;
118  f->m_fldToBnd = false;
119  }
120 
122  {
123  }
124 
125  void ProcessDisplacement::Process(po::variables_map &vm)
126  {
127  if (m_f->m_verbose)
128  {
129  if(m_f->m_comm->GetRank() == 0)
130  {
131  cout << "ProcessDisplacement: Calculating displacement..."
132  << endl;
133  }
134  }
135 
136  string toFile = m_config["to"].as<string>();
137 
138  if (toFile == "")
139  {
140  cout << "ProcessDisplacement: you must provide a file" << endl;
141  return;
142  }
143 
144  bool useVertexIds = m_config["usevertexids"].m_beenSet;
145 
146  vector<string> files;
147  files.push_back(toFile);
152 
153  // Try to find boundary condition expansion.
154  int bndCondId = m_config["id"].as<int>();
155 
156  // FIXME: We should be storing boundary condition IDs
157  // somewhere...
158  m_f->m_bndRegionsToWrite.push_back(bndCondId);
159 
160  if (bndGraph->GetMeshDimension() == 1)
161  {
162  m_f->m_exp.push_back(m_f->AppendExpList(0, "v"));
163 
164  MultiRegions::ExpListSharedPtr bndCondExpU =
165  m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
166  MultiRegions::ExpListSharedPtr bndCondExpV =
167  m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
168 
169  map<int, int> bndCondIds;
170  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
171  {
172  bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()]
173  = i;
174  }
175 
176  const SpatialDomains::SegGeomMap &tmp =
177  bndGraph->GetAllSegGeoms();
178  SpatialDomains::SegGeomMap::const_iterator sIt;
179 
180  for (sIt = tmp.begin(); sIt != tmp.end(); ++sIt)
181  {
182  map<int, int>::iterator mIt = bndCondIds.find(sIt->first);
183 
184  if (mIt == bndCondIds.end())
185  {
186  cout << "Warning: couldn't find element "
187  << sIt->first << endl;
188  continue;
189  }
190 
191  int e = mIt->second;
192 
194  boost::dynamic_pointer_cast<SpatialDomains::SegGeom>(
195  bndCondExpU->GetExp(e)->GetGeom());
196 
197  SpatialDomains::SegGeomSharedPtr to = sIt->second;
198 
199  // Create temporary SegExp
201  LocalRegions::SegExp>::AllocateSharedPtr(
202  bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
203  to);
204 
205  const int offset = bndCondExpU->GetPhys_Offset(e);
206  const int nq = toSeg->GetTotPoints();
207 
208  Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
209 
210  bndCondExpU->GetExp(e)->GetCoords(xC, yC);
211  toSeg->GetCoords(xL, yL);
212 
213  Vmath::Vsub(nq, xL, 1, xC, 1,
214  tmp = bndCondExpU->UpdatePhys() + offset, 1);
215  Vmath::Vsub(nq, yL, 1, yC, 1,
216  tmp = bndCondExpV->UpdatePhys() + offset, 1);
217  }
218 
219  // bndconstrained?
220  bndCondExpU->FwdTrans_BndConstrained(
221  bndCondExpU->GetPhys(), bndCondExpU->UpdateCoeffs());
222  bndCondExpV->FwdTrans_BndConstrained(
223  bndCondExpV->GetPhys(), bndCondExpV->UpdateCoeffs());
224  }
225  else if (bndGraph->GetMeshDimension() == 2)
226  {
227  m_f->m_exp.push_back(m_f->AppendExpList(0, "v"));
228  m_f->m_exp.push_back(m_f->AppendExpList(0, "w"));
229 
230  MultiRegions::ExpListSharedPtr bndCondExpU =
231  m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
232  MultiRegions::ExpListSharedPtr bndCondExpV =
233  m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
234  MultiRegions::ExpListSharedPtr bndCondExpW =
235  m_f->m_exp[2]->GetBndCondExpansions()[bndCondId];
236 
237  map<int, int> bndCondIds;
238  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
239  {
240  bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()]
241  = i;
242  }
243 
244  TriFaceMap vertexFaceMap;
245 
246  if (useVertexIds)
247  {
248  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
249  {
251  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
252  bndCondExpU->GetExp(i)->GetGeom());
253 
254  TriFaceIDs t(from->GetVid(0), from->GetVid(1),
255  from->GetVid(2));
256  vertexFaceMap[t] = i;
257  }
258  }
259 
260  const SpatialDomains::TriGeomMap &tmp =
261  bndGraph->GetAllTriGeoms();
262  SpatialDomains::TriGeomMap::const_iterator sIt;
263 
264  for (sIt = tmp.begin(); sIt != tmp.end(); ++sIt)
265  {
267  int e;
268 
269  if (useVertexIds)
270  {
271  TriFaceIDs t(sIt->second->GetVid(0),
272  sIt->second->GetVid(1),
273  sIt->second->GetVid(2));
274 
275  tIt = vertexFaceMap.find(t);
276  e = tIt == vertexFaceMap.end() ? -1 : tIt->second;
277  }
278  else
279  {
281  mIt = bndCondIds.find(sIt->first);
282  e = mIt == bndCondIds.end() ? -1 : mIt->second;
283  }
284 
285  if (e == -1)
286  {
287  cout << "Warning: couldn't find element "
288  << sIt->first << endl;
289  continue;
290  }
291 
293  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
294  bndCondExpU->GetExp(e)->GetGeom());
295 
296  SpatialDomains::TriGeomSharedPtr to = sIt->second;
297 
298  // Create temporary SegExp
300  LocalRegions::TriExp>::AllocateSharedPtr(
301  bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
302  bndCondExpV->GetExp(e)->GetBasis(1)->GetBasisKey(),
303  to);
304 
305  const int offset = bndCondExpU->GetPhys_Offset(e);
306  const int nq = toSeg->GetTotPoints();
307 
308  Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
309  Array<OneD, NekDouble> zL(nq), zC(nq);
310 
311  bndCondExpU->GetExp(e)->GetCoords(xC, yC, zC);
312  toSeg->GetCoords(xL, yL, zL);
313 
314  Vmath::Vsub(nq, xL, 1, xC, 1,
315  tmp = bndCondExpU->UpdatePhys() + offset, 1);
316  Vmath::Vsub(nq, yL, 1, yC, 1,
317  tmp = bndCondExpV->UpdatePhys() + offset, 1);
318  Vmath::Vsub(nq, zL, 1, zC, 1,
319  tmp = bndCondExpW->UpdatePhys() + offset, 1);
320  }
321 
322  // bndconstrained?
323  bndCondExpU->FwdTrans_BndConstrained(
324  bndCondExpU->GetPhys(), bndCondExpU->UpdateCoeffs());
325  bndCondExpV->FwdTrans_BndConstrained(
326  bndCondExpV->GetPhys(), bndCondExpV->UpdateCoeffs());
327  bndCondExpW->FwdTrans_BndConstrained(
328  bndCondExpW->GetPhys(), bndCondExpW->UpdateCoeffs());
329  }
330  }
331 }
332 }
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:121
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
virtual void Process()=0
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
FieldSharedPtr m_f
Field object.
std::size_t operator()(TriFaceIDs const &p) const
std::map< int, TriGeomSharedPtr > TriGeomMap
Definition: TriGeom.h:62
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
std::map< int, SegGeomSharedPtr > SegGeomMap
Definition: SegGeom.h:54
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
boost::unordered_map< TriFaceIDs, int, TriFaceHash > TriFaceMap
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.
Definition: Vmath.cpp:329
Represents a command-line configuration option.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
Definition: VtkToFld.cpp:53
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
boost::shared_ptr< TriExp > TriExpSharedPtr
Definition: TriExp.h:291
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215