BoundaryElementDist.cc 1.45 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include "BoundaryElementDist.h"

void 
BoundaryElementDist::calcNormal(
	             const FixVec<WorldVector<double>, DIMEN> &vecs, 
		     WorldVector<double> &normalVec)
{
  FUNCNAME("BoundaryElementDist::calcNormal");
  
  switch(dim) {
  case 2: calcNormal_2d(vecs, normalVec);
    break;
  case 3: calcNormal_3d(vecs, normalVec);
    break;
  default: ERROR_EXIT("illegal dimension !\n");
  }
}

void 
BoundaryElementDist::calcNormal_2d(
		     const FixVec<WorldVector<double>, DIMEN> &vecs, 
		     WorldVector<double> &normalVec)
{
  FUNCNAME("BoundaryElementDist::calcNormal_2d");
  
  TEST_EXIT(vecs.size() == dim)("illegal number of world vectors in vecs !\n");
  
  double norm2 = 0.0;
  double val;
  for (int i=0; i<dim; ++i){
    val = vecs[0][i] - vecs[1][i];
    norm2 += val*val;
    normalVec[dim-1-i] = val;
  } 
  normalVec[0] *= -1;
  double norm = sqrt(norm2);
  for(int i=0; i<dim; ++i) {
    normalVec[i] = 1/norm * normalVec[i];
  }
}

void 
BoundaryElementDist::calcNormal_3d(
		     const FixVec<WorldVector<double>, DIMEN> &vecs, 
		     WorldVector<double> &normalVec)
{
  FUNCNAME("BoundaryElementDist::calcNormal_3d");
  
  TEST_EXIT(vecs.size() == dim)("illegal number of world vectors in vecs !\n");
  
  WorldVector<double> A = vecs[1]-vecs[0];
  WorldVector<double> B = vecs[2]-vecs[0];
  vectorProduct(A, B, normalVec);
  
  double norm = sqrt(normalVec * normalVec);
  for(int i=0; i<dim; ++i) {
    normalVec[i] = 1/norm * normalVec[i];
  }
}