Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 <iostream>
37 #include <string>
38 using namespace std;
39 
40 #include "ProcessDisplacement.h"
41 
44 #include <LocalRegions/SegExp.h>
45 #include <LocalRegions/TriExp.h>
46 #include <StdRegions/StdQuadExp.h>
47 #include <StdRegions/StdSegExp.h>
48 #include <StdRegions/StdTriExp.h>
49 
50 #include <boost/tuple/tuple.hpp>
51 #include <boost/tuple/tuple_comparison.hpp>
52 #include <boost/unordered_set.hpp>
53 
54 namespace Nektar
55 {
56 namespace FieldUtils
57 {
58 struct TriFaceIDs
59 {
60  TriFaceIDs(int a, int b, int c) : a(a), b(b), c(c)
61  {
62  }
63  int a;
64  int b;
65  int c;
66 };
67 
68 struct TriFaceHash : std::unary_function<TriFaceIDs, std::size_t>
69 {
70  std::size_t operator()(TriFaceIDs const &p) const
71  {
72  std::vector<int> ids(3);
73 
74  ids[0] = p.a;
75  ids[1] = p.b;
76  ids[2] = p.c;
77 
78  std::sort(ids.begin(), ids.end());
79  return boost::hash_range(ids.begin(), ids.end());
80  }
81 };
82 
83 bool operator==(TriFaceIDs const &p1, TriFaceIDs const &p2)
84 {
85  std::vector<int> ids1(3), ids2(3);
86 
87  ids1[0] = p1.a;
88  ids1[1] = p1.b;
89  ids1[2] = p1.c;
90  ids2[0] = p2.a;
91  ids2[1] = p2.b;
92  ids2[2] = p2.c;
93 
94  std::sort(ids1.begin(), ids1.end());
95  std::sort(ids2.begin(), ids2.end());
96 
97  return ids1[0] == ids2[0] && ids1[1] == ids2[1] && ids1[2] == ids2[2];
98 }
99 
100 typedef boost::unordered_map<TriFaceIDs, int, TriFaceHash> TriFaceMap;
101 
102 ModuleKey ProcessDisplacement::className =
104  ModuleKey(eProcessModule, "displacement"),
105  ProcessDisplacement::create,
106  "Deform a mesh given an input field defining displacement");
107 
108 ProcessDisplacement::ProcessDisplacement(FieldSharedPtr f) : ProcessModule(f)
109 {
110  m_config["to"] =
111  ConfigOption(false, "", "Name of file containing high order boundary");
112  m_config["id"] =
113  ConfigOption(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->TreatAsRankZero())
130  {
131  cout << "ProcessDisplacement: Calculating displacement..." << endl;
132  }
133  }
134 
135  string toFile = m_config["to"].as<string>();
136 
137  if (toFile == "")
138  {
139  cout << "ProcessDisplacement: you must provide a file" << endl;
140  return;
141  }
142 
143  bool useVertexIds = m_config["usevertexids"].m_beenSet;
144 
145  vector<string> files;
146  files.push_back(toFile);
151 
152  // Try to find boundary condition expansion.
153  int bndCondId = m_config["id"].as<int>();
154 
155  // FIXME: We should be storing boundary condition IDs
156  // somewhere...
157  m_f->m_bndRegionsToWrite.push_back(bndCondId);
158 
159  if (bndGraph->GetMeshDimension() == 1)
160  {
161  m_f->m_exp.push_back(m_f->AppendExpList(0, "v"));
162 
163  MultiRegions::ExpListSharedPtr bndCondExpU =
164  m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
165  MultiRegions::ExpListSharedPtr bndCondExpV =
166  m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
167 
168  map<int, int> bndCondIds;
169  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
170  {
171  bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()] = i;
172  }
173 
174  const SpatialDomains::SegGeomMap &tmp = 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 " << sIt->first << endl;
184  continue;
185  }
186 
187  int e = mIt->second;
188 
190  boost::dynamic_pointer_cast<SpatialDomains::SegGeom>(
191  bndCondExpU->GetExp(e)->GetGeom());
192 
193  SpatialDomains::SegGeomSharedPtr to = sIt->second;
194 
195  // Create temporary SegExp
198  bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(), to);
199 
200  const int offset = bndCondExpU->GetPhys_Offset(e);
201  const int nq = toSeg->GetTotPoints();
202 
203  Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
204 
205  bndCondExpU->GetExp(e)->GetCoords(xC, yC);
206  toSeg->GetCoords(xL, yL);
207 
208  Vmath::Vsub(nq, xL, 1, xC, 1,
209  tmp = bndCondExpU->UpdatePhys() + offset, 1);
210  Vmath::Vsub(nq, yL, 1, yC, 1,
211  tmp = bndCondExpV->UpdatePhys() + offset, 1);
212  }
213 
214  // bndconstrained?
215  bndCondExpU->FwdTrans_BndConstrained(bndCondExpU->GetPhys(),
216  bndCondExpU->UpdateCoeffs());
217  bndCondExpV->FwdTrans_BndConstrained(bndCondExpV->GetPhys(),
218  bndCondExpV->UpdateCoeffs());
219  }
220  else if (bndGraph->GetMeshDimension() == 2)
221  {
222  m_f->m_exp.push_back(m_f->AppendExpList(0, "v"));
223  m_f->m_exp.push_back(m_f->AppendExpList(0, "w"));
224 
225  MultiRegions::ExpListSharedPtr bndCondExpU =
226  m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
227  MultiRegions::ExpListSharedPtr bndCondExpV =
228  m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
229  MultiRegions::ExpListSharedPtr bndCondExpW =
230  m_f->m_exp[2]->GetBndCondExpansions()[bndCondId];
231 
232  map<int, int> bndCondIds;
233  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
234  {
235  bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()] = i;
236  }
237 
238  TriFaceMap vertexFaceMap;
239 
240  if (useVertexIds)
241  {
242  for (int i = 0; i < bndCondExpU->GetExpSize(); ++i)
243  {
245  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
246  bndCondExpU->GetExp(i)->GetGeom());
247 
248  TriFaceIDs t(from->GetVid(0), from->GetVid(1), from->GetVid(2));
249  vertexFaceMap[t] = i;
250  }
251  }
252 
253  const SpatialDomains::TriGeomMap &tmp = bndGraph->GetAllTriGeoms();
254  SpatialDomains::TriGeomMap::const_iterator sIt;
255 
256  for (sIt = tmp.begin(); sIt != tmp.end(); ++sIt)
257  {
259  int e;
260 
261  if (useVertexIds)
262  {
263  TriFaceIDs t(sIt->second->GetVid(0), sIt->second->GetVid(1),
264  sIt->second->GetVid(2));
265 
266  tIt = vertexFaceMap.find(t);
267  e = tIt == vertexFaceMap.end() ? -1 : tIt->second;
268  }
269  else
270  {
272  mIt = bndCondIds.find(sIt->first);
273  e = mIt == bndCondIds.end() ? -1 : mIt->second;
274  }
275 
276  if (e == -1)
277  {
278  cout << "Warning: couldn't find element " << sIt->first << endl;
279  continue;
280  }
281 
283  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
284  bndCondExpU->GetExp(e)->GetGeom());
285 
286  SpatialDomains::TriGeomSharedPtr to = sIt->second;
287 
288  // Create temporary SegExp
291  bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
292  bndCondExpV->GetExp(e)->GetBasis(1)->GetBasisKey(), to);
293 
294  const int offset = bndCondExpU->GetPhys_Offset(e);
295  const int nq = toSeg->GetTotPoints();
296 
297  Array<OneD, NekDouble> xL(nq), xC(nq), yL(nq), yC(nq), tmp;
298  Array<OneD, NekDouble> zL(nq), zC(nq);
299 
300  bndCondExpU->GetExp(e)->GetCoords(xC, yC, zC);
301  toSeg->GetCoords(xL, yL, zL);
302 
303  Vmath::Vsub(nq, xL, 1, xC, 1,
304  tmp = bndCondExpU->UpdatePhys() + offset, 1);
305  Vmath::Vsub(nq, yL, 1, yC, 1,
306  tmp = bndCondExpV->UpdatePhys() + offset, 1);
307  Vmath::Vsub(nq, zL, 1, zC, 1,
308  tmp = bndCondExpW->UpdatePhys() + offset, 1);
309  }
310 
311  // bndconstrained?
312  bndCondExpU->FwdTrans_BndConstrained(bndCondExpU->GetPhys(),
313  bndCondExpU->UpdateCoeffs());
314  bndCondExpV->FwdTrans_BndConstrained(bndCondExpV->GetPhys(),
315  bndCondExpV->UpdateCoeffs());
316  bndCondExpW->FwdTrans_BndConstrained(bndCondExpW->GetPhys(),
317  bndCondExpW->UpdateCoeffs());
318  }
319 }
320 }
321 }
virtual void Process(po::variables_map &vm)
Write mesh to output file.
map< string, ConfigOption > m_config
List of configuration values.
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:124
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Represents a command-line configuration option.
STL namespace.
pair< ModuleType, string > ModuleKey
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
std::map< int, TriGeomSharedPtr > TriGeomMap
Definition: TriGeom.h:62
boost::unordered_map< TriFaceIDs, int, TriFaceHash > TriFaceMap
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:270
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:767
std::map< int, SegGeomSharedPtr > SegGeomMap
Definition: SegGeom.h:54
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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:343
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
Abstract base class for processing modules.
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:54
std::size_t operator()(TriFaceIDs const &p) const
boost::shared_ptr< TriExp > TriExpSharedPtr
Definition: TriExp.h:295
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215