50 #include <boost/tuple/tuple.hpp>
51 #include <boost/tuple/tuple_comparison.hpp>
52 #include <boost/unordered_set.hpp>
68 struct TriFaceHash : std::unary_function<TriFaceIDs, std::size_t>
72 std::vector<int> ids(3);
78 std::sort(ids.begin(), ids.end());
79 return boost::hash_range(ids.begin(), ids.end());
85 std::vector<int> ids1(3), ids2(3);
94 std::sort(ids1.begin(), ids1.end());
95 std::sort(ids2.begin(), ids2.end());
97 return ids1[0] == ids2[0] && ids1[1] == ids2[1] && ids1[2] == ids2[2];
100 typedef boost::unordered_map<TriFaceIDs, int, TriFaceHash>
TriFaceMap;
102 ModuleKey ProcessDisplacement::className =
105 ProcessDisplacement::create,
106 "Deform a mesh given an input field defining displacement");
111 ConfigOption(
false,
"",
"Name of file containing high order boundary");
113 ConfigOption(
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 if (
m_f->m_comm->TreatAsRankZero())
131 cout <<
"ProcessDisplacement: Calculating displacement..." << endl;
135 string toFile =
m_config[
"to"].as<
string>();
139 cout <<
"ProcessDisplacement: you must provide a file" << endl;
143 bool useVertexIds =
m_config[
"usevertexids"].m_beenSet;
145 vector<string> files;
146 files.push_back(toFile);
153 int bndCondId =
m_config[
"id"].as<
int>();
157 m_f->m_bndRegionsToWrite.push_back(bndCondId);
159 if (bndGraph->GetMeshDimension() == 1)
161 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"v"));
164 m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
166 m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
168 map<int, int> bndCondIds;
169 for (
int i = 0; i < bndCondExpU->GetExpSize(); ++i)
171 bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()] = i;
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 " << sIt->first << endl;
191 bndCondExpU->GetExp(e)->GetGeom());
198 bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(), to);
200 const int offset = bndCondExpU->GetPhys_Offset(e);
201 const int nq = toSeg->GetTotPoints();
205 bndCondExpU->GetExp(e)->GetCoords(xC, yC);
206 toSeg->GetCoords(xL, yL);
209 tmp = bndCondExpU->UpdatePhys() + offset, 1);
211 tmp = bndCondExpV->UpdatePhys() + offset, 1);
215 bndCondExpU->FwdTrans_BndConstrained(bndCondExpU->GetPhys(),
216 bndCondExpU->UpdateCoeffs());
217 bndCondExpV->FwdTrans_BndConstrained(bndCondExpV->GetPhys(),
218 bndCondExpV->UpdateCoeffs());
220 else if (bndGraph->GetMeshDimension() == 2)
222 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"v"));
223 m_f->m_exp.push_back(
m_f->AppendExpList(0,
"w"));
226 m_f->m_exp[0]->GetBndCondExpansions()[bndCondId];
228 m_f->m_exp[1]->GetBndCondExpansions()[bndCondId];
230 m_f->m_exp[2]->GetBndCondExpansions()[bndCondId];
232 map<int, int> bndCondIds;
233 for (
int i = 0; i < bndCondExpU->GetExpSize(); ++i)
235 bndCondIds[bndCondExpU->GetExp(i)->GetGeom()->GetGlobalID()] = i;
238 TriFaceMap vertexFaceMap;
242 for (
int i = 0; i < bndCondExpU->GetExpSize(); ++i)
246 bndCondExpU->GetExp(i)->GetGeom());
248 TriFaceIDs t(from->GetVid(0), from->GetVid(1), from->GetVid(2));
249 vertexFaceMap[t] = i;
254 SpatialDomains::TriGeomMap::const_iterator sIt;
256 for (sIt = tmp.begin(); sIt != tmp.end(); ++sIt)
263 TriFaceIDs t(sIt->second->GetVid(0), sIt->second->GetVid(1),
264 sIt->second->GetVid(2));
266 tIt = vertexFaceMap.find(t);
267 e = tIt == vertexFaceMap.end() ? -1 : tIt->second;
272 mIt = bndCondIds.find(sIt->first);
273 e = mIt == bndCondIds.end() ? -1 : mIt->second;
278 cout <<
"Warning: couldn't find element " << sIt->first << endl;
284 bndCondExpU->GetExp(e)->GetGeom());
291 bndCondExpU->GetExp(e)->GetBasis(0)->GetBasisKey(),
292 bndCondExpV->GetExp(e)->GetBasis(1)->GetBasisKey(), to);
294 const int offset = bndCondExpU->GetPhys_Offset(e);
295 const int nq = toSeg->GetTotPoints();
300 bndCondExpU->GetExp(e)->GetCoords(xC, yC, zC);
301 toSeg->GetCoords(xL, yL, zL);
304 tmp = bndCondExpU->UpdatePhys() + offset, 1);
306 tmp = bndCondExpV->UpdatePhys() + offset, 1);
308 tmp = bndCondExpW->UpdatePhys() + offset, 1);
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());
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)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Represents a command-line configuration option.
TriFaceIDs(int a, int b, int c)
pair< ModuleType, string > ModuleKey
virtual ~ProcessDisplacement()
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
std::map< int, TriGeomSharedPtr > TriGeomMap
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
boost::shared_ptr< SegGeom > SegGeomSharedPtr
boost::shared_ptr< Field > FieldSharedPtr
std::map< int, SegGeomSharedPtr > SegGeomMap
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.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Abstract base class for processing modules.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
std::size_t operator()(TriFaceIDs const &p) const
boost::shared_ptr< TriExp > TriExpSharedPtr
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.