Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_STK_SurfaceNodeNormals.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
44 
45 #include "Panzer_STK_Interface.hpp"
49 #include "Panzer_Workset_Builder.hpp"
52 #include "Panzer_CellData.hpp"
54 
55 #include <stk_mesh/base/Selector.hpp>
56 #include <stk_mesh/base/GetEntities.hpp>
57 #include <stk_mesh/base/GetBuckets.hpp>
58 #include <stk_mesh/base/CreateAdjacentEntities.hpp>
59 
60 #include "Shards_CellTopology.hpp"
61 //#include "Intrepid2_FunctionSpaceTools.hpp"
62 #include "Intrepid2_CellTools_Serial.hpp"
63 #include "Teuchos_Assert.hpp"
64 
65 namespace panzer_stk {
66 
67  void computeSidesetNodeNormals(std::unordered_map<unsigned,std::vector<double> >& normals,
69  const std::string& sidesetName,
70  const std::string& elementBlockName,
71  std::ostream* /* out */,
72  std::ostream* pout)
73  {
74  using panzer::Cell;
75  using panzer::NODE;
76  using panzer::Dim;
77 
78  using Teuchos::RCP;
79 
80  panzer::MDFieldArrayFactory af("",true);
81 
82  RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
83  RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
84 
85  // Grab all nodes for a surface including ghosted to get correct contributions to normal average
86  stk::mesh::Part * sidePart = mesh->getSideset(sidesetName);
87  stk::mesh::Part * elmtPart = mesh->getElementBlockPart(elementBlockName);
88  stk::mesh::Selector sideSelector = *sidePart;
89  stk::mesh::Selector blockSelector = *elmtPart;
90  stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
91  std::vector<stk::mesh::Entity> sides;
92  stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
93 
94  std::vector<std::size_t> localSideTopoIDs;
95  std::vector<stk::mesh::Entity> parentElements;
96  panzer_stk::workset_utils::getUniversalSubcellElements(*mesh,elementBlockName,sides,localSideTopoIDs,parentElements);
97 
98  if (pout != NULL) {
99  for (std::size_t i=0; i < localSideTopoIDs.size(); ++i) {
100  *pout << "parent element: "
101  << " gid(" << bulkData->identifier(parentElements[i]) << ")"
102  << ", local_face(" << localSideTopoIDs[i] << ")"
103  << std::endl;
104  }
105  }
106 
107  // Do a single element at a time so that we don't have to batch up
108  // into faces
109 
110  // maps a panzer local element id to a list of normals
111  std::unordered_map<unsigned,std::vector<double> > nodeNormals;
112 
113  TEUCHOS_ASSERT(sides.size() == localSideTopoIDs.size());
114  TEUCHOS_ASSERT(localSideTopoIDs.size() == parentElements.size());
115 
116  RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
117  //Intrepid2::DefaultCubatureFactory cubFactory;
118  int cubDegree = 1;
119 
120  std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
121  std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
122  std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
123 
124  // KK: invoke serial interface; cubDegree is 1 and integration point is one
125  // for debugging statement, use max dimension
126  // this lookup table setup is necessary before any impl::celltools is called
127  Intrepid2::Impl::CellTools::setSubcellParametrization();
128  Kokkos::DynRankView<double,Kokkos::HostSpace> normal_at_point("normal",3); // parentTopology->getDimension());
129  for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
130 
131  std::vector<stk::mesh::Entity> elementEntities;
132  elementEntities.push_back(*parentElement); // notice this is size 1!
133  PHX::MDField<double,panzer::Cell,panzer::NODE,panzer::Dim> vertices
134  = af.buildStaticArray<double,Cell,NODE,Dim>("",elementEntities.size(), parentTopology->getVertexCount(), mesh->getDimension());
135  mesh->getElementVerticesNoResize(elementEntities,elementBlockName,vertices);
136 
137  panzer::CellData sideCellData(1,*sideID,parentTopology); // this is size 1 because elementEntties is size 1!
138  RCP<panzer::IntegrationRule> ir = Teuchos::rcp(new panzer::IntegrationRule(cubDegree,sideCellData));
139 
141  iv.setupArrays(ir);
142  iv.evaluateValues(vertices);
143 
144  // KK: use serial interface; jac_at_point (D,D) from (C,P,D,D)
145  {
146  auto jac_at_point = Kokkos::subview(iv.jac.get_view(), 0, 0, Kokkos::ALL(), Kokkos::ALL());
147  Intrepid2::Impl::
148  CellTools::Serial::getPhysicalSideNormal(normal_at_point, jac_at_point, *sideID, *(ir->topology));
149  }
150  // Kokkos::DynRankView<double,PHX::Device> normal("normal",1,ir->num_points,parentTopology->getDimension());
151  // Intrepid2::CellTools<PHX::exec_space>::getPhysicalSideNormals(normal, iv.jac.get_view(), *sideID, *(ir->topology));
152 
153  if (pout != NULL) {
154  *pout << "element normals: "
155  << "gid(" << bulkData->identifier(*parentElement) << ")"
156  << ", normal(" << normal_at_point(0) << "," << normal_at_point(1) << "," << normal_at_point(2) << ")"
157  << std::endl;
158  }
159 
160  // loop over nodes in nodes in side and add normal contribution for averaging
161  const size_t numNodes = bulkData->num_nodes(*side);
162  stk::mesh::Entity const* nodeRelations = bulkData->begin_nodes(*side);
163  for (size_t n=0; n<numNodes; ++n) {
164  stk::mesh::Entity node = nodeRelations[n];
165  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
166  nodeNormals[bulkData->identifier(node)].push_back(normal_at_point(dim));
167  }
168  }
169 
170  }
171 
172  // Now do the averaging of contributions
173  //std::unordered_map<unsigned,std::vector<double> > normals;
174  for (std::unordered_map<unsigned,std::vector<double> >::const_iterator node = nodeNormals.begin(); node != nodeNormals.end(); ++node) {
175 
176  TEUCHOS_ASSERT( (node->second.size() % parentTopology->getDimension()) == 0);
177 
178  int numSurfaceContributions = node->second.size() / parentTopology->getDimension();
179  std::vector<double> contribArea(numSurfaceContributions);
180  double totalArea = 0.0;
181  for (int surface = 0; surface < numSurfaceContributions; ++surface) {
182 
183  double sum = 0.0;
184  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
185  sum += (node->second[surface*parentTopology->getDimension() + dim]) *
186  (node->second[surface*parentTopology->getDimension() + dim]);
187 
188  contribArea[surface] = std::sqrt(sum);
189 
190  totalArea += contribArea[surface];
191  }
192 
193  // change the contribArea to the scale factor for each contribution
194  for (std::size_t i = 0; i < contribArea.size(); ++i)
195  contribArea[i] /= totalArea;
196 
197  // loop over contributions and compute final normal values
198  normals[node->first].resize(parentTopology->getDimension());
199  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
200  normals[node->first][dim] = 0.0;
201  for (int surface = 0; surface < numSurfaceContributions; ++surface) {
202  normals[node->first][dim] += node->second[surface*parentTopology->getDimension() + dim] * contribArea[surface] / totalArea;
203  }
204  }
205 
206  if (pout != NULL) {
207  *pout << "surface normal before normalization: "
208  << "gid(" << node->first << ")"
209  << ", normal(" << normals[node->first][0] << "," << normals[node->first][1] << "," << normals[node->first][2] << ")"
210  << std::endl;
211  }
212 
213  double sum = 0.0;
214  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
215  sum += normals[node->first][dim] * normals[node->first][dim];
216 
217  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
218  normals[node->first][dim] /= std::sqrt(sum);
219 
220  if (pout != NULL) {
221  *pout << "surface normal after normalization: "
222  << "gid(" << node->first << ")"
223  << ", normal("
224  << normals[node->first][0] << ","
225  << normals[node->first][1] << ","
226  << normals[node->first][2] << ")"
227  << std::endl;
228  }
229 
230  }
231 
232  }
233 
234  void computeSidesetNodeNormals(std::unordered_map<std::size_t,Kokkos::DynRankView<double,PHX::Device> >& normals,
236  const std::string& sidesetName,
237  const std::string& elementBlockName,
238  std::ostream* out,
239  std::ostream* pout)
240  {
241  using Teuchos::RCP;
242 
243  std::unordered_map<unsigned,std::vector<double> > nodeEntityIdToNormals;
244 
245  computeSidesetNodeNormals(nodeEntityIdToNormals,mesh,sidesetName,elementBlockName,out,pout);
246 
247  RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
248  RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
249 
250  // Grab all nodes for a surface including ghosted to get correct contributions to normal average
251  stk::mesh::Part * sidePart = mesh->getSideset(sidesetName);
252  stk::mesh::Part * elmtPart = mesh->getElementBlockPart(elementBlockName);
253  stk::mesh::Selector sideSelector = *sidePart;
254  stk::mesh::Selector blockSelector = *elmtPart;
255  stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
256  std::vector<stk::mesh::Entity> sides;
257  stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
258 
259  RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
260 
261  std::vector<std::size_t> localSideTopoIDs;
262  std::vector<stk::mesh::Entity> parentElements;
263  panzer_stk::workset_utils::getUniversalSubcellElements(*mesh,elementBlockName,sides,localSideTopoIDs,parentElements);
264 
265  std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
266  std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
267  std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
268  for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
269 
270  // loop over nodes in nodes in side element
271  const size_t numNodes = bulkData->num_nodes(*parentElement);
272  stk::mesh::Entity const* nodeRelations = bulkData->begin_nodes(*parentElement);
273 
274  normals[mesh->elementLocalId(*parentElement)] = Kokkos::DynRankView<double,PHX::Device>("normals",numNodes,parentTopology->getDimension());
275 
276  for (size_t nodeIndex=0; nodeIndex<numNodes; ++nodeIndex) {
277  stk::mesh::Entity node = nodeRelations[nodeIndex];
278  // if the node is on the sideset, insert, otherwise set normal
279  // to zero (it is an interior node of the parent element).
280  if (nodeEntityIdToNormals.find(bulkData->identifier(node)) != nodeEntityIdToNormals.end()) {
281  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
282  (normals[mesh->elementLocalId(*parentElement)])(nodeIndex,dim) = (nodeEntityIdToNormals[bulkData->identifier(node)])[dim];
283  }
284  }
285  else {
286  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
287  (normals[mesh->elementLocalId(*parentElement)])(nodeIndex,dim) = 0.0;
288  }
289  }
290  }
291 
292  }
293 
294  }
295 
296 }
void computeSidesetNodeNormals(std::unordered_map< unsigned, std::vector< double > > &normals, const Teuchos::RCP< const panzer_stk::STK_Interface > &mesh, const std::string &sidesetName, const std::string &elementBlockName, std::ostream *, std::ostream *pout)
Computes the normals for all nodes associated with a sideset surface.
void getElementVerticesNoResize(const std::vector< std::size_t > &localIds, ArrayT &vertices) const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block count
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
stk::mesh::Part * getSideset(const std::string &name) const
std::size_t elementLocalId(stk::mesh::Entity elmt) const
unsigned getDimension() const
get the dimension
void getUniversalSubcellElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &entities, std::vector< std::size_t > &localEntityIds, std::vector< stk::mesh::Entity > &elements)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Data for determining cell topology and dimensionality.
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int num_cells=-1)
Evaluate basis values.
Teuchos::RCP< const shards::CellTopology > topology
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
BASIS NODE
Spatial Dimension Tag.
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const