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 // Panzer: A partial differential equation assembly
4 // engine for strongly coupled complex multiphysics systems
5 //
6 // Copyright 2011 NTESS and the Panzer contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
12 
13 #include "Panzer_STK_Interface.hpp"
17 #include "Panzer_Workset_Builder.hpp"
20 #include "Panzer_CellData.hpp"
22 
23 #include <stk_mesh/base/Selector.hpp>
24 #include <stk_mesh/base/GetEntities.hpp>
25 #include <stk_mesh/base/CreateAdjacentEntities.hpp>
26 
27 #include "Shards_CellTopology.hpp"
28 //#include "Intrepid2_FunctionSpaceTools.hpp"
29 #include "Intrepid2_CellTools_Serial.hpp"
30 #include "Teuchos_Assert.hpp"
31 
32 namespace panzer_stk {
33 
34  void computeSidesetNodeNormals(std::unordered_map<unsigned,std::vector<double> >& normals,
36  const std::string& sidesetName,
37  const std::string& elementBlockName,
38  std::ostream* /* out */,
39  std::ostream* pout)
40  {
41  using panzer::Cell;
42  using panzer::NODE;
43  using panzer::Dim;
44 
45  using Teuchos::RCP;
46 
47  panzer::MDFieldArrayFactory af("",true);
48 
49  RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
50  RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
51 
52  // Grab all nodes for a surface including ghosted to get correct contributions to normal average
53  stk::mesh::Part * sidePart = mesh->getSideset(sidesetName);
54  stk::mesh::Part * elmtPart = mesh->getElementBlockPart(elementBlockName);
55  stk::mesh::Selector sideSelector = *sidePart;
56  stk::mesh::Selector blockSelector = *elmtPart;
57  stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
58  std::vector<stk::mesh::Entity> sides;
59  stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
60 
61  std::vector<std::size_t> missingElementIndices;
62  std::vector<std::size_t> localSideTopoIDs;
63  std::vector<stk::mesh::Entity> parentElements;
64  panzer_stk::workset_utils::getUniversalSubcellElements(*mesh,elementBlockName,sides,localSideTopoIDs,parentElements,missingElementIndices);
65 
66  // Delete all sides whose neighboring element in elementBlockName is not in the current process
67  std::vector<std::size_t>::reverse_iterator index;
68  for(index=missingElementIndices.rbegin(); index != missingElementIndices.rend(); ++index) {
69  sides.erase(sides.begin()+*index);
70  }
71 
72  if (pout != NULL) {
73  for (std::size_t i=0; i < localSideTopoIDs.size(); ++i) {
74  *pout << "parent element: "
75  << " gid(" << bulkData->identifier(parentElements[i]) << ")"
76  << ", local_face(" << localSideTopoIDs[i] << ")"
77  << std::endl;
78  }
79  }
80 
81  // Do a single element at a time so that we don't have to batch up
82  // into faces
83 
84  // maps a panzer local element id to a list of normals
85  std::unordered_map<unsigned,std::vector<double> > nodeNormals;
86 
87  TEUCHOS_ASSERT(sides.size() == localSideTopoIDs.size());
88  TEUCHOS_ASSERT(localSideTopoIDs.size() == parentElements.size());
89 
90  RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
91  //Intrepid2::DefaultCubatureFactory cubFactory;
92  int cubDegree = 1;
93 
94  std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
95  std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
96  std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
97 
98  // KK: invoke serial interface; cubDegree is 1 and integration point is one
99  // for debugging statement, use max dimension
100  auto side_parametrization = Intrepid2::RefSubcellParametrization<Kokkos::HostSpace>::get(2,parentTopology->getKey());
101  Kokkos::DynRankView<double,Kokkos::HostSpace> normal_at_point("normal",3); // parentTopology->getDimension());
102  for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
103 
104  std::vector<stk::mesh::Entity> elementEntities;
105  elementEntities.push_back(*parentElement); // notice this is size 1!
107  = af.buildStaticArray<double,Cell,NODE,Dim>("",elementEntities.size(), parentTopology->getNodeCount(), mesh->getDimension());
108  auto node_view = nodes.get_view();
109  mesh->getElementNodesNoResize(elementEntities,elementBlockName,node_view);
110 
111  panzer::CellData sideCellData(1,*sideID,parentTopology); // this is size 1 because elementEntties is size 1!
112  RCP<panzer::IntegrationRule> ir = Teuchos::rcp(new panzer::IntegrationRule(cubDegree,sideCellData));
113 
115  iv.setupArrays(ir);
116  iv.evaluateValues(nodes);
117 
118  // KK: use serial interface; jac_at_point (D,D) from (C,P,D,D)
119  {
120  auto jac_at_point = Kokkos::subview(iv.jac.get_view(), 0, 0, Kokkos::ALL(), Kokkos::ALL());
121  auto jac_at_point_h = Kokkos::create_mirror_view(jac_at_point);
122  Kokkos::deep_copy(jac_at_point_h, jac_at_point);
123  Intrepid2::Impl::
124  CellTools::Serial::getPhysicalSideNormal(normal_at_point, side_parametrization, jac_at_point_h, *sideID);
125  }
126 
127  if (pout != NULL) {
128  *pout << "element normals: "
129  << "gid(" << bulkData->identifier(*parentElement) << ")"
130  << ", normal(" << normal_at_point(0) << "," << normal_at_point(1) << "," << normal_at_point(2) << ")"
131  << std::endl;
132  }
133 
134  // loop over nodes in nodes in side and add normal contribution for averaging
135  const size_t numNodes = bulkData->num_nodes(*side);
136  stk::mesh::Entity const* nodeRelations = bulkData->begin_nodes(*side);
137  for (size_t n=0; n<numNodes; ++n) {
138  stk::mesh::Entity node = nodeRelations[n];
139  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
140  nodeNormals[bulkData->identifier(node)].push_back(normal_at_point(dim));
141  }
142  }
143 
144  }
145 
146  // Now do the averaging of contributions
147  //std::unordered_map<unsigned,std::vector<double> > normals;
148  for (std::unordered_map<unsigned,std::vector<double> >::const_iterator node = nodeNormals.begin(); node != nodeNormals.end(); ++node) {
149 
150  TEUCHOS_ASSERT( (node->second.size() % parentTopology->getDimension()) == 0);
151 
152  int numSurfaceContributions = node->second.size() / parentTopology->getDimension();
153  std::vector<double> contribArea(numSurfaceContributions);
154  double totalArea = 0.0;
155  for (int surface = 0; surface < numSurfaceContributions; ++surface) {
156 
157  double sum = 0.0;
158  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
159  sum += (node->second[surface*parentTopology->getDimension() + dim]) *
160  (node->second[surface*parentTopology->getDimension() + dim]);
161 
162  contribArea[surface] = std::sqrt(sum);
163 
164  totalArea += contribArea[surface];
165  }
166 
167  // change the contribArea to the scale factor for each contribution
168  for (std::size_t i = 0; i < contribArea.size(); ++i)
169  contribArea[i] /= totalArea;
170 
171  // loop over contributions and compute final normal values
172  normals[node->first].resize(parentTopology->getDimension());
173  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
174  normals[node->first][dim] = 0.0;
175  for (int surface = 0; surface < numSurfaceContributions; ++surface) {
176  normals[node->first][dim] += node->second[surface*parentTopology->getDimension() + dim] * contribArea[surface] / totalArea;
177  }
178  }
179 
180  if (pout != NULL) {
181  *pout << "surface normal before normalization: "
182  << "gid(" << node->first << ")"
183  << ", normal(" << normals[node->first][0] << "," << normals[node->first][1] << "," << normals[node->first][2] << ")"
184  << std::endl;
185  }
186 
187  double sum = 0.0;
188  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
189  sum += normals[node->first][dim] * normals[node->first][dim];
190 
191  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
192  normals[node->first][dim] /= std::sqrt(sum);
193 
194  if (pout != NULL) {
195  *pout << "surface normal after normalization: "
196  << "gid(" << node->first << ")"
197  << ", normal("
198  << normals[node->first][0] << ","
199  << normals[node->first][1] << ","
200  << normals[node->first][2] << ")"
201  << std::endl;
202  }
203 
204  }
205 
206  }
207 
208  void computeSidesetNodeNormals(std::unordered_map<std::size_t,Kokkos::DynRankView<double,PHX::Device> >& normals,
210  const std::string& sidesetName,
211  const std::string& elementBlockName,
212  std::ostream* out,
213  std::ostream* pout)
214  {
215  using Teuchos::RCP;
216 
217  std::unordered_map<unsigned,std::vector<double> > nodeEntityIdToNormals;
218 
219  computeSidesetNodeNormals(nodeEntityIdToNormals,mesh,sidesetName,elementBlockName,out,pout);
220 
221  RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
222  RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
223 
224  // Grab all nodes for a surface including ghosted to get correct contributions to normal average
225  stk::mesh::Part * sidePart = mesh->getSideset(sidesetName);
226  stk::mesh::Part * elmtPart = mesh->getElementBlockPart(elementBlockName);
227  stk::mesh::Selector sideSelector = *sidePart;
228  stk::mesh::Selector blockSelector = *elmtPart;
229  stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
230  std::vector<stk::mesh::Entity> sides;
231  stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
232 
233  RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
234 
235  std::vector<std::size_t> missingElementIndices;
236  std::vector<std::size_t> localSideTopoIDs;
237  std::vector<stk::mesh::Entity> parentElements;
238  panzer_stk::workset_utils::getUniversalSubcellElements(*mesh,elementBlockName,sides,localSideTopoIDs,parentElements,missingElementIndices);
239 
240  // Delete all sides whose neighboring element in elementBlockName is not in the current process
241  std::vector<std::size_t>::reverse_iterator index;
242  for(index=missingElementIndices.rbegin(); index != missingElementIndices.rend(); ++index) {
243  sides.erase(sides.begin()+*index);
244  }
245 
246  std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
247  std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
248  std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
249  for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
250 
251  // loop over nodes in nodes in side element
252  const size_t numNodes = bulkData->num_nodes(*parentElement);
253  stk::mesh::Entity const* nodeRelations = bulkData->begin_nodes(*parentElement);
254 
255  normals[mesh->elementLocalId(*parentElement)] = Kokkos::DynRankView<double,PHX::Device>("normals",numNodes,parentTopology->getDimension());
256  auto normals_h = Kokkos::create_mirror_view(normals[mesh->elementLocalId(*parentElement)]);
257  for (size_t nodeIndex=0; nodeIndex<numNodes; ++nodeIndex) {
258  stk::mesh::Entity node = nodeRelations[nodeIndex];
259  // if the node is on the sideset, insert, otherwise set normal
260  // to zero (it is an interior node of the parent element).
261  if (nodeEntityIdToNormals.find(bulkData->identifier(node)) != nodeEntityIdToNormals.end()) {
262  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
263  normals_h(nodeIndex,dim) = (nodeEntityIdToNormals[bulkData->identifier(node)])[dim];
264  }
265  }
266  else {
267  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
268  normals_h(nodeIndex,dim) = 0.0;
269  }
270  }
271  }
272  Kokkos::deep_copy(normals[mesh->elementLocalId(*parentElement)], normals_h);
273  }
274 
275  }
276 
277 }
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 getElementNodesNoResize(const std::vector< std::size_t > &localIds, ArrayT &nodes) const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
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
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Data for determining cell topology and dimensionality.
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, std::vector< std::size_t > &missingElementIndices)
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
#define TEUCHOS_ASSERT(assertion_test)
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &cell_node_coordinates, const int num_cells=-1, const Teuchos::RCP< const SubcellConnectivity > &face_connectivity=Teuchos::null, const int num_virtual_cells=-1)
Evaluate basis values.
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const