Intrepid2
Intrepid2_TensorTopologyMap.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid2 Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov),
38// Mauro Perego (mperego@sandia.gov), or
39// Nate Roberts (nvrober@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49#ifndef Intrepid2_TensorTopologyMap_h
50#define Intrepid2_TensorTopologyMap_h
51
52#include "Intrepid2_Utils.hpp"
53
54#include "Shards_CellTopology.hpp"
55
56#include <map>
57#include <set>
58#include <vector>
59
60namespace Intrepid2
61{
71 {
72 shards::CellTopology cellTopo1_;
73 shards::CellTopology cellTopo2_;
74 shards::CellTopology compositeCellTopo_;
75
76 using Subcell = std::pair<unsigned,unsigned>; // first: dimension, second: subcell ordinal
77 using SubcellPair = std::pair<Subcell,Subcell>; // first: subcell in cellTopo1_, second: subcell in cellTopo2_
78 std::map<SubcellPair,Subcell> subcellMap_; // values are the subcell in compositeCellTopo (the dimension of this subcell is the sum of the dimensions in the SubcellPair)
79 public:
94 TensorTopologyMap(shards::CellTopology cellTopo1, shards::CellTopology cellTopo2, shards::CellTopology compositeCellTopo,
95 std::vector< std::pair<unsigned,unsigned> > nodePairs = std::vector< std::pair<unsigned,unsigned> >())
96 :
97 cellTopo1_(cellTopo1),
98 cellTopo2_(cellTopo2),
99 compositeCellTopo_(compositeCellTopo)
100 {
101 if (nodePairs.size() == 0)
102 {
103 nodePairs = defaultNodePairs(cellTopo1, cellTopo2, compositeCellTopo);
104 }
105
106 auto spaceDim1 = cellTopo1.getDimension();
107 auto spaceDim2 = cellTopo2.getDimension();
108 auto spaceDim = compositeCellTopo.getDimension();
109 INTREPID2_TEST_FOR_EXCEPTION(spaceDim1 + spaceDim2 != spaceDim, std::invalid_argument, "incompatible spatial dimensions");
110 std::map<std::pair<unsigned,unsigned>,unsigned> compositeNodeOrdinalMap;
111 for (unsigned compositeNodeOrdinal=0; compositeNodeOrdinal<nodePairs.size(); compositeNodeOrdinal++)
112 {
113 compositeNodeOrdinalMap[nodePairs[compositeNodeOrdinal]] = compositeNodeOrdinal;
114 }
115 for (unsigned d1=0; d1<=spaceDim1; d1++)
116 {
117 unsigned subcellCount1 = cellTopo1.getSubcellCount(d1);
118 for (unsigned subcellOrdinal1=0; subcellOrdinal1<subcellCount1; subcellOrdinal1++)
119 {
120 Subcell subcell1 = {d1, subcellOrdinal1};
121 std::set<unsigned> subcell1Nodes; // set because we don't care about ordering
122 unsigned nodeCount1 = cellTopo1_.getNodeCount(d1, subcellOrdinal1);
123 // unfortunately, the node count for vertices is given as 0 by shards::CellTopology. This seems like a bug.
124 if (d1 == 0)
125 {
126 subcell1Nodes.insert(subcellOrdinal1);
127 }
128 else
129 {
130 for (unsigned nodeOrdinal1=0; nodeOrdinal1<nodeCount1; nodeOrdinal1++)
131 {
132 subcell1Nodes.insert(cellTopo1_.getNodeMap(d1, subcellOrdinal1, nodeOrdinal1));
133 }
134 }
135 for (unsigned d2=0; d2<=spaceDim2; d2++)
136 {
137 unsigned subcellCount2 = cellTopo2.getSubcellCount(d2);
138 for (unsigned subcellOrdinal2=0; subcellOrdinal2<subcellCount2; subcellOrdinal2++)
139 {
140 Subcell subcell2 = {d2, subcellOrdinal2};
141 std::set<unsigned> subcell2Nodes; // set because we don't care about ordering
142 unsigned nodeCount2 = cellTopo2_.getNodeCount(d2, subcellOrdinal2);
143 // unfortunately, the node count for vertices is given as 0 by shards::CellTopology. This seems like a bug.
144 if (d2 == 0)
145 {
146 subcell2Nodes.insert(subcellOrdinal2);
147 }
148 else
149 {
150 for (unsigned nodeOrdinal2=0; nodeOrdinal2<nodeCount2; nodeOrdinal2++)
151 {
152 subcell2Nodes.insert(cellTopo2_.getNodeMap(d2, subcellOrdinal2, nodeOrdinal2));
153 }
154 }
155
156 std::set<unsigned> compositeNodes; // all the nodes from subcell1 times nodes from subcell2
157 for (auto subcellNode1 : subcell1Nodes)
158 {
159 for (auto subcellNode2 : subcell2Nodes)
160 {
161 INTREPID2_TEST_FOR_EXCEPTION(compositeNodeOrdinalMap.find({subcellNode1,subcellNode2}) == compositeNodeOrdinalMap.end(),
162 std::invalid_argument, "Node combination not found in map");
163 compositeNodes.insert(compositeNodeOrdinalMap[{subcellNode1,subcellNode2}]);
164 }
165 }
166 // now, search the composite topology for the unique subcell that involves all those nodes
167 // we do a brute-force search; we'll never have very big topologies, so the cost is not an issue
168 unsigned compositeSubcellDim = d1 + d2;
169 unsigned compositeSubcellCount = compositeCellTopo.getSubcellCount(compositeSubcellDim);
170 bool compositeSubcellFound = false;
171 for (unsigned compositeSubcellOrdinal=0; compositeSubcellOrdinal<compositeSubcellCount; compositeSubcellOrdinal++)
172 {
173 unsigned compositeSubcellNodeCount = (compositeSubcellDim > 0) ? compositeCellTopo.getNodeCount(compositeSubcellDim, compositeSubcellOrdinal)
174 : 1; // again, dealing with the fact that the node count for vertices is defined as 0, not 1
175 if (compositeSubcellNodeCount != compositeNodes.size()) continue; // node counts don't match, so this is not the subcell we're looking for
176 bool matches = true; // if we don't find a node that isn't contained in compositeNodes, then this is the subcell we're looking for
177 for (unsigned compositeSubcellNode=0; compositeSubcellNode<compositeSubcellNodeCount; compositeSubcellNode++)
178 {
179 unsigned nodeInCell = compositeCellTopo.getNodeMap(compositeSubcellDim, compositeSubcellOrdinal, compositeSubcellNode);
180 if (compositeNodes.find(nodeInCell) == compositeNodes.end())
181 {
182 matches = false;
183 break;
184 }
185 }
186 if (matches)
187 {
188 compositeSubcellFound = true;
189 subcellMap_[{subcell1,subcell2}] = {compositeSubcellDim, compositeSubcellOrdinal};
190 break;
191 }
192 }
193 INTREPID2_TEST_FOR_EXCEPTION(!compositeSubcellFound, std::invalid_argument, "Composite subcell not found");
194 }
195 }
196 }
197 }
198 }
199 TensorTopologyMap(shards::CellTopology cellTopo1, shards::CellTopology cellTopo2)
200 :
201 TensorTopologyMap(cellTopo1, cellTopo2, compositeCellTopology(cellTopo1,cellTopo2)) {}
202
214 unsigned getCompositeSubcellOrdinal(unsigned subcell1Dim, unsigned subcell1Ordinal, unsigned subcell2Dim, unsigned subcell2Ordinal)
215 {
216 Subcell subcell1 = {subcell1Dim, subcell1Ordinal};
217 Subcell subcell2 = {subcell2Dim, subcell2Ordinal};
218 auto entry = subcellMap_.find({subcell1,subcell2});
219 INTREPID2_TEST_FOR_EXCEPTION(entry == subcellMap_.end(), std::invalid_argument, "entry not found");
220 auto subcell = entry->second;
221 return subcell.second; // subcell ordinal
222 }
223
236 static std::vector< std::vector<int> > referenceCellGeometry(shards::CellTopology cellTopo)
237 {
238 std::vector< std::vector<int> > nodes(cellTopo.getVertexCount());
239 auto key = cellTopo.getKey();
240 switch (key)
241 {
242 case shards::Node::key:
243 nodes.push_back({}); // node.getVertexCount() returns 0; we do want a single (empty) entry, though, even though there's no spatial variation
244 break;
245 case shards::Line<2>::key:
246 nodes[0] = {0};
247 nodes[1] = {1};
248 break;
249 case shards::Triangle<3>::key:
250 nodes[0] = {0,0};
251 nodes[1] = {1,0};
252 nodes[2] = {0,1};
253 break;
254 case shards::Quadrilateral<4>::key:
255 nodes[0] = {0,0};
256 nodes[1] = {1,0};
257 nodes[2] = {1,1};
258 nodes[3] = {0,1};
259 break;
260 case shards::Hexahedron<8>::key:
261 nodes[0] = {0,0,0};
262 nodes[1] = {1,0,0};
263 nodes[2] = {1,1,0};
264 nodes[3] = {0,1,0};
265 nodes[4] = {0,0,1};
266 nodes[5] = {1,0,1};
267 nodes[6] = {1,1,1};
268 nodes[7] = {0,1,1};
269 break;
270 case shards::Wedge<6>::key: // wedge is triangle in (x,y) times line in z
271 nodes[0] = {0,0,0};
272 nodes[1] = {1,0,0};
273 nodes[2] = {0,1,0};
274 nodes[3] = {0,0,1};
275 nodes[4] = {1,0,1};
276 nodes[5] = {0,1,1};
277 break;
278 default:
279 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported CellTopology");
280 }
281 return nodes;
282 }
283
297 static std::vector< std::pair<unsigned,unsigned> > defaultNodePairs(shards::CellTopology cellTopo1, shards::CellTopology cellTopo2, shards::CellTopology compositeCellTopo)
298 {
299 unsigned compositeNodeCount = compositeCellTopo.getVertexCount();
300 unsigned nodeCount1 = cellTopo1.getVertexCount();
301 unsigned nodeCount2 = cellTopo2.getVertexCount();
302 INTREPID2_TEST_FOR_EXCEPTION(compositeNodeCount != nodeCount1 * nodeCount2, std::invalid_argument, "Incompatible topologies");
303 std::vector< std::pair<unsigned,unsigned> > nodePairs(compositeNodeCount);
304 auto nodes1 = referenceCellGeometry(cellTopo1);
305 auto nodes2 = referenceCellGeometry(cellTopo2);
306 auto compositeNodes = referenceCellGeometry(compositeCellTopo);
307 std::map< std::vector<int>, unsigned > compositeNodeMap;
308 for (unsigned compositeNodeOrdinal=0; compositeNodeOrdinal<compositeNodeCount; compositeNodeOrdinal++)
309 {
310 compositeNodeMap[compositeNodes[compositeNodeOrdinal]] = compositeNodeOrdinal;
311 }
312 for (unsigned nodeOrdinal1=0; nodeOrdinal1<nodeCount1; nodeOrdinal1++)
313 {
314 std::vector<int> node1 = nodes1[nodeOrdinal1];
315 for (unsigned nodeOrdinal2=0; nodeOrdinal2<nodeCount2; nodeOrdinal2++)
316 {
317 std::vector<int> node2 = nodes2[nodeOrdinal2];
318 std::vector<int> compositeNode(node1); // copy node1 into the first slots of compositeNode
319 compositeNode.insert(compositeNode.end(), node2.begin(), node2.end());
320 auto compositeNodeMapEntry = compositeNodeMap.find(compositeNode);
321 INTREPID2_TEST_FOR_EXCEPTION(compositeNodeMapEntry == compositeNodeMap.end(), std::invalid_argument, "composite node not found in map");
322 nodePairs[compositeNodeMapEntry->second] = {nodeOrdinal1,nodeOrdinal2};
323 }
324 }
325 return nodePairs;
326 }
327
335 static shards::CellTopology compositeCellTopology(shards::CellTopology cellTopo1, shards::CellTopology cellTopo2)
336 {
337 if (cellTopo1.getBaseKey() == shards::Node::key)
338 {
339 return cellTopo2;
340 }
341 else if (cellTopo2.getBaseKey() == shards::Node::key)
342 {
343 return cellTopo1;
344 }
345
346 // if we get here, neither is a Node
347 auto key1 = cellTopo1.getBaseKey();
348 auto key2 = cellTopo2.getBaseKey();
349 switch (key1)
350 {
351 case shards::Line<2>::key:
352 switch (key2)
353 {
354 case shards::Line<2>::key:
355 return shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<> >() );
356 case shards::Triangle<3>::key:
357 // unsupported:
358 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"line x triangle is not supported at present. Wedge is triangle x line.");
359 case shards::Quadrilateral<4>::key:
360 return shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<> >() );
361 default:
362 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unsupported cell topology pair");
363 }
364 case shards::Triangle<3>::key:
365 switch (key2)
366 {
367 case shards::Line<2>::key:
368 return shards::CellTopology(shards::getCellTopologyData<shards::Wedge<> >() );
369 default:
370 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unsupported cell topology pair");
371 }
372 case shards::Quadrilateral<4>::key:
373 switch (key2)
374 {
375 case shards::Line<2>::key:
376 return shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<> >() );
377 default:
378 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unsupported cell topology pair");
379 }
380 default:
381 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unsupported cell topology pair");
382 }
383 }
384 };
385} // end namespace Intrepid2
386
387#endif /* Intrepid2_TensorTopology_h */
Header function for Intrepid2::Util class and other utility functions.
For two cell topologies whose tensor product is a third, this class establishes a mapping from subcel...
unsigned getCompositeSubcellOrdinal(unsigned subcell1Dim, unsigned subcell1Ordinal, unsigned subcell2Dim, unsigned subcell2Ordinal)
Map from component subcell ordinals to the corresponding composite subcell ordinal.
static std::vector< std::vector< int > > referenceCellGeometry(shards::CellTopology cellTopo)
A topological representation of the vertex geometries corresponding to a cell topology.
TensorTopologyMap(shards::CellTopology cellTopo1, shards::CellTopology cellTopo2, shards::CellTopology compositeCellTopo, std::vector< std::pair< unsigned, unsigned > > nodePairs=std::vector< std::pair< unsigned, unsigned > >())
Constructor.
static shards::CellTopology compositeCellTopology(shards::CellTopology cellTopo1, shards::CellTopology cellTopo2)
The composite cell topology generated as the tensor product of the specified component cell topologie...
static std::vector< std::pair< unsigned, unsigned > > defaultNodePairs(shards::CellTopology cellTopo1, shards::CellTopology cellTopo2, shards::CellTopology compositeCellTopo)
Default node pairs corresponding to the provided component and composite cell topologies.