Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessJac.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessJac.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Calculate Jacobians of elements.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 #include "ProcessJac.h"
38 
39 using namespace std;
40 using namespace Nektar::NekMeshUtils;
41 
42 namespace Nektar
43 {
44 namespace Utilities
45 {
46 
47 ModuleKey ProcessJac::className = GetModuleFactory().RegisterCreatorFunction(
48  ModuleKey(eProcessModule, "jac"),
49  ProcessJac::create,
50  "Process elements based on values of Jacobian.");
51 
52 ProcessJac::ProcessJac(MeshSharedPtr m) : ProcessModule(m)
53 {
54  m_config["extract"] =
55  ConfigOption(false, "0.0", "Extract non-valid elements from mesh.");
56  m_config["list"] = ConfigOption(
57  false, "0", "Print list of elements having negative Jacobian.");
58 }
59 
61 {
62 }
63 
65 {
66  if (m_mesh->m_verbose)
67  {
68  cout << "ProcessJac: Calculating Jacobians... " << endl;
69  }
70 
71  bool extract = m_config["extract"].beenSet;
72  bool printList = m_config["list"].as<bool>();
73  NekDouble thres = m_config["extract"].as<NekDouble>();
74 
75  vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
76 
77  if (extract)
78  {
79  m_mesh->m_element[m_mesh->m_expDim].clear();
80  }
81 
82  if (printList)
83  {
84  cout << "Elements with negative Jacobian:" << endl;
85  }
86 
87  int nNeg = 0;
88 
89  Array<OneD, int> bin(20, 0);
90 
91  // Iterate over list of elements of expansion dimension.
92  for (int i = 0; i < el.size(); ++i)
93  {
94  // Create elemental geometry.
96  el[i]->GetGeom(m_mesh->m_spaceDim);
97 
98  // Generate geometric factors.
99  SpatialDomains::GeomFactorsSharedPtr gfac = geom->GetGeomFactors();
100 
101  LibUtilities::PointsKeyVector p = geom->GetPointsKeys();
102  SpatialDomains::DerivStorage deriv = gfac->GetDeriv(p);
103  const int pts = deriv[0][0].num_elements();
104  Array<OneD,NekDouble> jc(pts);
105  for (int k = 0; k < pts; ++k)
106  {
107  DNekMat jac(m_mesh->m_expDim, m_mesh->m_expDim, 0.0, eFULL);
108 
109  for (int l = 0; l < m_mesh->m_expDim; ++l)
110  {
111  for (int j = 0; j < m_mesh->m_expDim; ++j)
112  {
113  jac(j,l) = deriv[l][j][k];
114  }
115  }
116 
117  if(m_mesh->m_expDim == 2)
118  {
119  jc[k] = jac(0,0) * jac(1,1) - jac(0,1)*jac(1,0);
120  }
121  else if(m_mesh->m_expDim == 3)
122  {
123  jc[k] = jac(0,0) * (jac(1,1)*jac(2,2) - jac(2,1)*jac(1,2)) -
124  jac(0,1) * (jac(1,0)*jac(2,2) - jac(2,0)*jac(1,2)) +
125  jac(0,2) * (jac(1,0)*jac(2,1) - jac(2,0)*jac(1,1));
126  }
127  }
128 
129  NekDouble scaledJac = Vmath::Vmin(jc.num_elements(),jc,1) /
130  Vmath::Vmax(jc.num_elements(),jc,1);
131 
132  bool valid = gfac->IsValid();
133 
134  if (extract && (scaledJac < thres || !valid))
135  {
136  m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
137  }
138 
139  // Get the Jacobian and, if it is negative, print a warning
140  // message.
141  if (!valid)
142  {
143  nNeg++;
144 
145  if (printList)
146  {
147  cout << " - " << el[i]->GetId() << " ("
148  << LibUtilities::ShapeTypeMap[el[i]->GetConf().m_e] << ")"
149  << " " << scaledJac
150  << endl;
151  }
152  }
153  }
154 
155  if (extract)
156  {
157  m_mesh->m_element[m_mesh->m_expDim - 1].clear();
158  ProcessVertices();
159  ProcessEdges();
160  ProcessFaces();
161  ProcessElements();
163  }
164 
165  if (printList || m_mesh->m_verbose)
166  {
167  cout << "Total negative Jacobians: " << nNeg << endl;
168  }
169  else if (nNeg > 0)
170  {
171  cout << "WARNING: Detected " << nNeg << " element"
172  << (nNeg == 1 ? "" : "s") << " with negative Jacobian." << endl;
173  }
174 }
175 }
176 }
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
virtual void Process()
Write mesh to output file.
Definition: ProcessJac.cpp:64
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:779
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:871
STL namespace.
pair< ModuleType, string > ModuleKey
const char *const ShapeTypeMap[]
Definition: ShapeType.hpp:66
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
Represents a command-line configuration option.
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
Abstract base class for processing modules.
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
std::pair< ModuleType, std::string > ModuleKey
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215