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  cout << "ProcessDisplacement: Calculating displacement..."
130  << endl;
131  }
132 
133  string toFile = m_config["to"].as<string>();
134 
135  if (toFile == "")
136  {
137  cout << "ProcessDisplacement: you must provide a file" << endl;
138  return;
139  }
140 
141  bool useVertexIds = m_config["usevertexids"].m_beenSet;
142 
143  vector<string> files;
144  files.push_back(toFile);
149 
150  // Try to find boundary condition expansion.
151  int bndCondId = m_config["id"].as<int>();
152 
153  // FIXME: We should be storing boundary condition IDs
154  // somewhere...
155  m_f->m_bndRegionsToWrite.push_back(bndCondId);
156 
157  if (bndGraph->GetMeshDimension() == 1)
158  {
159  m_f->m_exp.push_back(m_f->AppendExpList(0, "v"));
160 
161  MultiRegions::ExpListSharedPtr bndCondExpU =
162  m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
163  MultiRegions::ExpListSharedPtr bndCondExpV =
164  m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
165 
166  map<int, int> bndCondIds;
167  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
168  {
169  bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()]
170  = i;
171  }
172 
173  const SpatialDomains::SegGeomMap &tmp =
174  bndGraph->GetAllSegGeoms();
175  SpatialDomains::SegGeomMap::const_iterator sIt;
176 
177  for (sIt = tmp.begin(); sIt != tmp.end(); ++sIt)
178  {
179  map<int, int>::iterator mIt = bndCondIds.find(sIt->first);
180 
181  if (mIt == bndCondIds.end())
182  {
183  cout << "Warning: couldn't find element "
184  << sIt->first << endl;
185  continue;
186  }
187 
188  int e = mIt->second;
189 
191  boost::dynamic_pointer_cast<SpatialDomains::SegGeom>(
192  bndCondExpU->GetExp(e)->GetGeom());
193 
194  SpatialDomains::SegGeomSharedPtr to = sIt->second;
195 
196  // Create temporary SegExp
198  LocalRegions::SegExp>::AllocateSharedPtr(
199  bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
200  to);
201 
202  const int offset = bndCondExpU->GetPhys_Offset(e);
203  const int nq = toSeg->GetTotPoints();
204 
205  Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
206 
207  bndCondExpU->GetExp(e)->GetCoords(xC, yC);
208  toSeg->GetCoords(xL, yL);
209 
210  Vmath::Vsub(nq, xL, 1, xC, 1,
211  tmp = bndCondExpU->UpdatePhys() + offset, 1);
212  Vmath::Vsub(nq, yL, 1, yC, 1,
213  tmp = bndCondExpV->UpdatePhys() + offset, 1);
214  }
215 
216  // bndconstrained?
217  bndCondExpU->FwdTrans_BndConstrained(
218  bndCondExpU->GetPhys(), bndCondExpU->UpdateCoeffs());
219  bndCondExpV->FwdTrans_BndConstrained(
220  bndCondExpV->GetPhys(), bndCondExpV->UpdateCoeffs());
221  }
222  else if (bndGraph->GetMeshDimension() == 2)
223  {
224  m_f->m_exp.push_back(m_f->AppendExpList(0, "v"));
225  m_f->m_exp.push_back(m_f->AppendExpList(0, "w"));
226 
227  MultiRegions::ExpListSharedPtr bndCondExpU =
228  m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
229  MultiRegions::ExpListSharedPtr bndCondExpV =
230  m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
231  MultiRegions::ExpListSharedPtr bndCondExpW =
232  m_f->m_exp[2]->GetBndCondExpansions()[bndCondId];
233 
234  map<int, int> bndCondIds;
235  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
236  {
237  bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()]
238  = i;
239  }
240 
241  TriFaceMap vertexFaceMap;
242 
243  if (useVertexIds)
244  {
245  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
246  {
248  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
249  bndCondExpU->GetExp(i)->GetGeom());
250 
251  TriFaceIDs t(from->GetVid(0), from->GetVid(1),
252  from->GetVid(2));
253  vertexFaceMap[t] = i;
254  }
255  }
256 
257  const SpatialDomains::TriGeomMap &tmp =
258  bndGraph->GetAllTriGeoms();
259  SpatialDomains::TriGeomMap::const_iterator sIt;
260 
261  for (sIt = tmp.begin(); sIt != tmp.end(); ++sIt)
262  {
264  int e;
265 
266  if (useVertexIds)
267  {
268  TriFaceIDs t(sIt->second->GetVid(0),
269  sIt->second->GetVid(1),
270  sIt->second->GetVid(2));
271 
272  tIt = vertexFaceMap.find(t);
273  e = tIt == vertexFaceMap.end() ? -1 : tIt->second;
274  }
275  else
276  {
278  mIt = bndCondIds.find(sIt->first);
279  e = mIt == bndCondIds.end() ? -1 : mIt->second;
280  }
281 
282  if (e == -1)
283  {
284  cout << "Warning: couldn't find element "
285  << sIt->first << endl;
286  continue;
287  }
288 
290  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
291  bndCondExpU->GetExp(e)->GetGeom());
292 
293  SpatialDomains::TriGeomSharedPtr to = sIt->second;
294 
295  // Create temporary SegExp
297  LocalRegions::TriExp>::AllocateSharedPtr(
298  bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
299  bndCondExpV->GetExp(e)->GetBasis(1)->GetBasisKey(),
300  to);
301 
302  const int offset = bndCondExpU->GetPhys_Offset(e);
303  const int nq = toSeg->GetTotPoints();
304 
305  Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
306  Array<OneD, NekDouble> zL(nq), zC(nq);
307 
308  bndCondExpU->GetExp(e)->GetCoords(xC, yC, zC);
309  toSeg->GetCoords(xL, yL, zL);
310 
311  Vmath::Vsub(nq, xL, 1, xC, 1,
312  tmp = bndCondExpU->UpdatePhys() + offset, 1);
313  Vmath::Vsub(nq, yL, 1, yC, 1,
314  tmp = bndCondExpV->UpdatePhys() + offset, 1);
315  Vmath::Vsub(nq, zL, 1, zC, 1,
316  tmp = bndCondExpW->UpdatePhys() + offset, 1);
317  }
318 
319  // bndconstrained?
320  bndCondExpU->FwdTrans_BndConstrained(
321  bndCondExpU->GetPhys(), bndCondExpU->UpdateCoeffs());
322  bndCondExpV->FwdTrans_BndConstrained(
323  bndCondExpV->GetPhys(), bndCondExpV->UpdateCoeffs());
324  bndCondExpW->FwdTrans_BndConstrained(
325  bndCondExpW->GetPhys(), bndCondExpW->UpdateCoeffs());
326  }
327  }
328 }
329 }
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