math_lib.cc 3.87 KB
Newer Older
Phillip Lino Rall's avatar
Phillip Lino Rall committed
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
/**********************************************************************************
 * Copyright 2010 Christoph Pflaum 
 * 		Department Informatik Lehrstuhl 10 - Systemsimulation
 *		Friedrich-Alexander Universität Erlangen-Nürnberg
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 **********************************************************************************/




#include "../abbrevi.h"
#include "../parameter.h"
#include "math_lib.h"

// ------------------------------------------------------------
// math_lib.cc
//
// ------------------------------------------------------------

////////////////////////////////////
//   Find next prime 
////////////////////////////////////

double calc_maximal_edge_angle(const D3vector& va,
			       const D3vector& vb,
			       const D3vector& vc,
			       const D3vector& vd) {
  return MAX(MAX(max_interior_angel_of_triangle(vb,vc,vd),
		 max_interior_angel_of_triangle(va,vc,vb)),
	     MAX(max_interior_angel_of_triangle(va,vd,vc),
		 max_interior_angel_of_triangle(va,vb,vd)));
}

double calc_maximal_face_angle(const D3vector& v0,
			       const D3vector& v1,
			       const D3vector& v2,
			       const D3vector& v3) {
  return MAX(MAX(angle_between_faces(v1,v0,v2,v3),
		 angle_between_faces(v2,v0,v3,v1),
		 angle_between_faces(v3,v0,v1,v2)),
	     MAX(angle_between_faces(v0,v2,v1,v3),
		 angle_between_faces(v0,v3,v2,v1),
		 angle_between_faces(v0,v1,v3,v2)));
}

57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
D3matrix rotationMatrix(D3vector vIn, D3vector vOut)
{
  //calculates rotation matrix, which rotates vIn, such that is points into direction vOut
    if (D3VectorNorm(vIn) == 0 || D3VectorNorm(vOut) == 0)
    {
        return D3matrix();
    }
    D3vector cross = cross_product(vOut,vIn);
    double normed = D3VectorNorm(cross_product(vOut,vIn));
    D3vector orthogonal = cross / normed;
    double alpha = atan2(normed,product(vOut,vIn));

    double s = sin(alpha);
    double c = cos(alpha);
    double mc = 1.0 - c;
    double x = orthogonal.x;
    double y = orthogonal.y;
    double z  = orthogonal.z;


    D3vector cx(c + x * x * mc,x * y * mc + z * s,x * z * mc - y * s); // first column
    D3vector cy(x * y * mc - z * s,c + y * y * mc,y * z * mc + x * s); // second column
    D3vector cz(x * z * mc + y * s,y * z * mc - x * s,c + z * z * mc); // third column

    D3matrix ret(cx,cy,cz);
    ret.transpose();
    return ret;
}
Phillip Lino Rall's avatar
Phillip Lino Rall committed
85
86
87
88
89
90
91
92
93
94
95
96

////////////////////////////////////
//    lambda of p in tet
////////////////////////////////////


D3vector lambda_of_p_in_tet(D3vector p,
			    D3vector cA, D3vector cB,
			    D3vector cC, D3vector cD) {
  D3matrix M(cB-cA,cC-cA,cD-cA);

  D3vector t;
97
  //t = M.invert_apply(p-cA);
98
  //cout<<"calculating lambda = ";t.Print();cout<<endl;
Phillip Lino Rall's avatar
Phillip Lino Rall committed
99
100
101
102
  return M.invert_apply(p-cA);
}


103

Phillip Lino Rall's avatar
Phillip Lino Rall committed
104
105
106
107
////////////////////////////////////
//    D3vector
////////////////////////////////////

108
109


Phillip Lino Rall's avatar
Phillip Lino Rall committed
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
bool D3vector::testZwischen(D3vector A, D3vector B) { // testet obPunkt auf Gerade zwischen A und B ist.
	if(((A.x <= x && x <= B.x) ||
            (A.x >= x && x >= B.x))==false) {
	  cout << " x-test: Punkt nicht auf Gerade: ";
          cout << "\n A: "; A.Print();
          cout << "\n B: "; B.Print();
          cout << "\n v: "; Print();  
	  return false;
	}
	return true;
	/*
		D3matrix testAufGerade(v_test_A,vcurrent,v_test_B);
		double det = testAufGerade.Determinante();
		if(ABS(det) > 0.001) {
		 cout << " Punkt nicht auf Gerade: " << det << endl;  vcurrent.Print(); cout << endl; 
		}
		*/

}