← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to Mesh/qualityMeasures.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-02 18:12:15 UTC
  • mfrom: (1.2.8 upstream)
  • Revision ID: james.westby@ubuntu.com-20090902181215-yla8zvcas2ucvkm9
[Christophe Prud'homme]
* New upstream release
  + fixed surface mesh orientation bug introduced in 2.4.0;
  + mesh and graphics code refactoring;
  + small usability enhancements and bug fixes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Mesh/qualityMeasures.cpp
6
6
#include "qualityMeasures.h"
7
7
#include "BDS.h"
8
8
#include "MVertex.h"
9
 
#include "MElement.h"
 
9
#include "MTriangle.h"
 
10
#include "MTetrahedron.h"
10
11
#include "Numeric.h"
11
12
#include "FunctionSpace.h"
12
13
#include "GmshMessage.h"
66
67
    // condition number
67
68
  case QMTRI_COND:
68
69
    {
 
70
      /*
69
71
      double a [3] = {xc - xa, yc - ya, zc - za};
70
72
      double b [3] = {xb - xa, yb - ya, zb - za};
71
73
      double c [3] ; prodve(a, b, c); norme(c);
72
74
      double A[3][3] = {{a[0] , b[0] , c[0]} ,
73
75
                        {a[1] , b[1] , c[1]} ,
74
76
                        {a[2] , b[2] , c[2]}};
 
77
      */
75
78
      quality = -1;
76
79
    }
77
80
    break;
102
105
             const double &x4, const double &y4, const double &z4, 
103
106
             const gmshQualityMeasure4Tet &cr, double *volume)
104
107
{
105
 
  double quality;
106
108
  switch(cr){
107
109
  case QMTET_ONE:
108
110
    return 1.0;
173
175
  }
174
176
}
175
177
 
176
 
 
177
178
double mesh_functional_distorsion(MTriangle *t, double u, double v)
178
179
{
179
180
  // compute uncurved element jacobian d_u x and d_v x
180
181
  double mat[3][3];  
181
 
  double d1 = t->getPrimaryJacobian(u,v,0,mat);
182
 
  
183
 
   //double d1 = t->getJacobian(u,v,0,mat);
 
182
  t->getPrimaryJacobian(u, v, 0, mat);
 
183
  // t->getJacobian(u,v,0,mat);
184
184
  double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
185
185
  double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
186
186
  double normal1[3];
187
187
  prodve(v1, v2, normal1);
188
188
  double nn = sqrt(normal1[0]*normal1[0] + 
189
 
                   normal1[1]*normal1[1] + 
190
 
                   normal1[2]*normal1[2]);
 
189
                   normal1[1]*normal1[1] + 
 
190
                   normal1[2]*normal1[2]);
191
191
  
192
192
  // compute uncurved element jacobian d_u x and d_v x
193
193
  
194
 
  double d2 = t->getJacobian(u, v, 0, mat);
195
 
  
 
194
  t->getJacobian(u, v, 0, mat);
196
195
  double v1b[3] = {mat[0][0], mat[0][1], mat[0][2]};
197
196
  double v2b[3] = {mat[1][0], mat[1][1], mat[1][2]};
198
197
  double normal[3];
199
198
  prodve(v1b, v2b, normal);
200
199
  
201
 
  double sign;
 
200
  double sign = 1.0;
202
201
  prosca(normal1, normal, &sign);
203
202
  double det = norm3(normal) * (sign > 0 ? 1. : -1.) / nn;  
 
203
  //double det = norm3(normal);
204
204
 
205
205
  // compute distorsion
206
 
  double dist = std::min(1. / det, det); 
207
 
  return dist;
 
206
  //  double dist = std::min(1. / det, det); 
 
207
  return det;
208
208
}
209
209
 
 
210
double mesh_functional_distorsion_p2(MTriangle *t)
 
211
{
 
212
  double d = mesh_functional_distorsion(t,0.,0.);
 
213
  d = std::min(d,mesh_functional_distorsion(t,1.,0.));
 
214
  d = std::min(d,mesh_functional_distorsion(t,0.,1.));
 
215
  d = std::min(d,mesh_functional_distorsion(t,.5,0.));
 
216
  d = std::min(d,mesh_functional_distorsion(t,0.,0.5));
 
217
  d = std::min(d,mesh_functional_distorsion(t,.5,0.5));
 
218
  return d;
 
219
}
210
220
 
211
221
double qmDistorsionOfMapping (MTriangle *e)
212
222
{
213
223
  //  return 1.0;
214
 
  if (e->getPolynomialOrder() == 1)return 1.0;
 
224
  if (e->getPolynomialOrder() == 1) return 1.0;
 
225
  //  if (e->getPolynomialOrder() == 2) return mesh_functional_distorsion_p2(e);
 
226
 
215
227
  IntPt *pts;
216
228
  int npts;
217
229
  e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
218
230
  double dmin;
219
 
  for (int i=0;i<npts;i++){
 
231
  for (int i = 0 ; i < npts; i++){
220
232
    const double u = pts[i].pt[0];
221
233
    const double v = pts[i].pt[1];
222
 
    const double w = pts[i].pt[2];
223
 
    const double di  = mesh_functional_distorsion (e,u,v);
224
 
    dmin = (i==0)? di : std::min(dmin,di);
 
234
    const double di  = mesh_functional_distorsion (e, u, v);
 
235
    dmin = (i == 0)? di : std::min(dmin, di);
225
236
  }
226
237
  const gmshMatrix<double>& points = e->getFunctionSpace()->points;
227
238
 
228
 
  for (int i=0;i<e->getNumPrimaryVertices();i++) {
229
 
    const double u = points(i,0);
230
 
    const double v = points(i,1);
231
 
    const double di  = mesh_functional_distorsion (e,u,v);
232
 
    dmin = std::min(dmin,di);
 
239
  for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
 
240
    const double u = points(i, 0);
 
241
    const double v = points(i, 1);
 
242
    const double di  = mesh_functional_distorsion (e, u, v);
 
243
    dmin = std::min(dmin, di);
233
244
  }
234
245
  return dmin;
235
246
}
238
249
{
239
250
  // compute uncurved element jacobian d_u x and d_v x
240
251
  double mat[3][3];  
241
 
  t->getPrimaryJacobian(u,v,w, mat);
 
252
  t->getPrimaryJacobian(u, v, w, mat);
242
253
  
243
254
  const double det1 = det3x3(mat);
244
255
 
245
256
   //const double det1 = t->getJacobian(u,v,w,mat);
246
257
  // const double det1 = det3x3(mat);
247
 
  t->getJacobian(u,v,w,mat);
 
258
  t->getJacobian(u, v, w, mat);
248
259
  const double detN = det3x3(mat);
249
260
  // const double detN = det3x3(mat);
250
261
 
251
262
  //  printf("%g %g %g = %g %g\n",u,v,w,det1,detN);
252
263
 
253
264
  if (det1 == 0 || detN == 0) return 0;
254
 
  double dist = std::min(detN/det1, det1/detN); 
 
265
  double dist = std::min(detN / det1, det1 / detN); 
255
266
  return dist;
256
267
}
257
268
 
258
 
 
259
 
double qmDistorsionOfMapping (MTetrahedron *e)
 
269
double qmDistorsionOfMapping(MTetrahedron *e)
260
270
{
261
 
  if (e->getPolynomialOrder() == 1)return 1.0;
 
271
  if (e->getPolynomialOrder() == 1) return 1.0;
262
272
  IntPt *pts;
263
273
  int npts;
264
 
  e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
 
274
  e->getIntegrationPoints(e->getPolynomialOrder(), &npts, &pts);
265
275
  double dmin;
266
 
  for (int i=0;i<npts;i++){
 
276
  for (int i = 0; i < npts; i++){
267
277
    const double u = pts[i].pt[0];
268
278
    const double v = pts[i].pt[1];
269
279
    const double w = pts[i].pt[2];
270
 
    const double di  = mesh_functional_distorsion (e,u,v,w);
271
 
    dmin = (i==0)? di : std::min(dmin,di);
 
280
    const double di  = mesh_functional_distorsion(e, u, v, w);
 
281
    dmin = (i == 0) ? di : std::min(dmin, di);
272
282
  }
273
283
  
274
284
  const gmshMatrix<double>& points = e->getFunctionSpace()->points;
275
285
 
276
 
  for (int i=0;i<e->getNumPrimaryVertices();i++) {
277
 
    const double u = points(i,0);
278
 
    const double v = points(i,1);
279
 
    const double w = points(i,2);
280
 
    const double di  = mesh_functional_distorsion (e,u,v,w);
281
 
    dmin = std::min(dmin,di);
 
286
  for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
 
287
    const double u = points(i, 0);
 
288
    const double v = points(i, 1);
 
289
    const double w = points(i, 2);
 
290
    const double di  = mesh_functional_distorsion(e, u, v, w);
 
291
    dmin = std::min(dmin, di);
282
292
  }
283
 
  //  printf("DMIN = %g\n\n",dmin);
284
293
 
285
 
  return dmin< 0 ? 0 :dmin;
 
294
  return dmin;
286
295
}