50 #include <boost/unordered_set.hpp>
51 #include <boost/tuple/tuple.hpp>
52 #include <boost/tuple/tuple_comparison.hpp>
66 struct TriFaceHash : std::unary_function<TriFaceIDs, std::size_t>
70 std::vector<int> ids(3);
76 std::sort(ids.begin(), ids.end());
77 return boost::hash_range(ids.begin(), ids.end());
83 std::vector<int> ids1(3), ids2(3);
92 std::sort(ids1.begin(), ids1.end());
93 std::sort(ids2.begin(), ids2.end());
95 return ids1[0] == ids2[0] &&
100 typedef boost::unordered_map<TriFaceIDs, int, TriFaceHash>
TriFaceMap;
102 ModuleKey ProcessDisplacement::className =
105 "Deform a mesh given an input field defining displacement");
111 false,
"",
"Name of file containing high order boundary");
113 false,
"",
"Boundary ID to calculate displacement for");
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;
129 cout <<
"ProcessDisplacement: Calculating displacement..."
133 string toFile =
m_config[
"to"].as<
string>();
137 cout <<
"ProcessDisplacement: you must provide a file" << endl;
141 bool useVertexIds =
m_config[
"usevertexids"].m_beenSet;
143 vector<string> files;
144 files.push_back(toFile);
151 int bndCondId =
m_config[
"id"].as<
int>();
155 m_f->m_bndRegionsToWrite.push_back(bndCondId);
157 if (bndGraph->GetMeshDimension() == 1)
159 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"v"));
162 m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
164 m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
166 map<int, int> bndCondIds;
167 for (
int i = 0; i < bndCondExpU->GetExpSize(); ++i)
169 bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()]
174 bndGraph->GetAllSegGeoms();
175 SpatialDomains::SegGeomMap::const_iterator sIt;
177 for (sIt = tmp.begin(); sIt != tmp.end(); ++sIt)
181 if (mIt == bndCondIds.end())
183 cout <<
"Warning: couldn't find element "
184 << sIt->first << endl;
192 bndCondExpU->GetExp(e)->GetGeom());
199 bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
202 const int offset = bndCondExpU->GetPhys_Offset(e);
207 bndCondExpU->GetExp(e)->GetCoords(xC, yC);
208 toSeg->GetCoords(xL, yL);
211 tmp = bndCondExpU->UpdatePhys() + offset, 1);
213 tmp = bndCondExpV->UpdatePhys() + offset, 1);
217 bndCondExpU->FwdTrans_BndConstrained(
218 bndCondExpU->GetPhys(), bndCondExpU->UpdateCoeffs());
219 bndCondExpV->FwdTrans_BndConstrained(
220 bndCondExpV->GetPhys(), bndCondExpV->UpdateCoeffs());
222 else if (bndGraph->GetMeshDimension() == 2)
224 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"v"));
225 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"w"));
228 m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
230 m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
232 m_f->m_exp[2]->GetBndCondExpansions()[bndCondId];
234 map<int, int> bndCondIds;
235 for (
int i = 0; i < bndCondExpU->GetExpSize(); ++i)
237 bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()]
241 TriFaceMap vertexFaceMap;
245 for (
int i = 0; i < bndCondExpU->GetExpSize(); ++i)
249 bndCondExpU->GetExp(i)->GetGeom());
251 TriFaceIDs t(from->GetVid(0), from->GetVid(1),
253 vertexFaceMap[t] = i;
258 bndGraph->GetAllTriGeoms();
259 SpatialDomains::TriGeomMap::const_iterator sIt;
261 for (sIt = tmp.begin(); sIt != tmp.end(); ++sIt)
269 sIt->second->GetVid(1),
270 sIt->second->GetVid(2));
272 tIt = vertexFaceMap.find(t);
273 e = tIt == vertexFaceMap.end() ? -1 : tIt->second;
278 mIt = bndCondIds.find(sIt->first);
279 e = mIt == bndCondIds.end() ? -1 : mIt->second;
284 cout <<
"Warning: couldn't find element "
285 << sIt->first << endl;
291 bndCondExpU->GetExp(e)->GetGeom());
298 bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
299 bndCondExpV->GetExp(e)->GetBasis(1)->GetBasisKey(),
302 const int offset = bndCondExpU->GetPhys_Offset(e);
308 bndCondExpU->GetExp(e)->GetCoords(xC, yC, zC);
309 toSeg->GetCoords(xL, yL, zL);
312 tmp = bndCondExpU->UpdatePhys() + offset, 1);
314 tmp = bndCondExpV->UpdatePhys() + offset, 1);
316 tmp = bndCondExpW->UpdatePhys() + offset, 1);
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());
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
map< string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
FieldSharedPtr m_f
Field object.
std::size_t operator()(TriFaceIDs const &p) const
std::map< int, TriGeomSharedPtr > TriGeomMap
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
boost::shared_ptr< SegExp > SegExpSharedPtr
boost::shared_ptr< SegGeom > SegGeomSharedPtr
std::map< int, SegGeomSharedPtr > SegGeomMap
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
boost::shared_ptr< Field > FieldSharedPtr
boost::unordered_map< TriFaceIDs, int, TriFaceHash > TriFaceMap
virtual ~ProcessDisplacement()
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.
Represents a command-line configuration option.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
TriFaceIDs(int a, int b, int c)
boost::shared_ptr< TriGeom > TriGeomSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
boost::shared_ptr< TriExp > TriExpSharedPtr
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.