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
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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
// Copyright (c) 2010-2026, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-443271.
//
// This file is part of the GLVis visualization tool and library. For more
// information and source code availability see https://glvis.org.
//
// GLVis is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifndef GLVIS_GEOM_UTILS_HPP
#define GLVIS_GEOM_UTILS_HPP
#include <mfem.hpp>
// Some inline functions
inline void LinearCombination(const double a, const double x[],
const double b, const double y[], double z[])
{
z[0] = a*x[0] + b*y[0];
z[1] = a*x[1] + b*y[1];
z[2] = a*x[2] + b*y[2];
}
inline double InnerProd(const double a[], const double b[])
{
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
inline void CrossProd(const double a[], const double b[], double cp[])
{
cp[0] = a[1] * b[2] - a[2] * b[1];
cp[1] = a[2] * b[0] - a[0] * b[2];
cp[2] = a[0] * b[1] - a[1] * b[0];
}
inline int Normalize(double v[])
{
double len = sqrt(InnerProd(v, v));
if (len > 0.0)
{
len = 1.0 / len;
}
else
{
return 1;
}
for (int i = 0; i < 3; i++)
{
v[i] *= len;
}
return 0;
}
inline int Normalize(mfem::DenseMatrix &normals)
{
int err = 0;
for (int i = 0; i < normals.Width(); i++)
{
err += Normalize(&normals(0, i));
}
return err;
}
// 'v' must define three vectors that add up to the zero vector
// that is v[0], v[1], and v[2] are the sides of a triangle
inline int UnitCrossProd(double v[][3], double nor[])
{
// normalize the three vectors
for (int i = 0; i < 3; i++)
if (Normalize(v[i]))
{
return 1;
}
// find the pair that forms an angle closest to pi/2, i.e. having
// the longest cross product:
double cp[3][3], max_a = 0.;
int k = 0;
for (int i = 0; i < 3; i++)
{
CrossProd(v[(i+1)%3], v[(i+2)%3], cp[i]);
double a = sqrt(InnerProd(cp[i], cp[i]));
if (max_a < a)
{
max_a = a, k = i;
}
}
if (max_a == 0.)
{
return 1;
}
for (int i = 0; i < 3; i++)
{
nor[i] = cp[k][i] / max_a;
}
return 0;
}
inline int Compute3DUnitNormal(const double p1[], const double p2[],
const double p3[], double nor[])
{
double v[3][3];
for (int i = 0; i < 3; i++)
{
v[0][i] = p2[i] - p1[i];
v[1][i] = p3[i] - p2[i];
v[2][i] = p1[i] - p3[i];
}
return UnitCrossProd(v, nor);
}
inline int Compute3DUnitNormal (const double p1[], const double p2[],
const double p3[], const double p4[], double nor[])
{
double v[3][3];
for (int i = 0; i < 3; i++)
{
// cross product of the two diagonals:
/*
v[0][i] = p3[i] - p1[i];
v[1][i] = p4[i] - p2[i];
v[2][i] = (p1[i] + p2[i]) - (p3[i] + p4[i]);
*/
// cross product of the two vectors connecting the midpoints of the
// two pairs of opposing sides; this gives the normal vector in the
// midpoint of the quad:
v[0][i] = 0.5 * ((p2[i] + p3[i]) - (p1[i] + p4[i]));
v[1][i] = 0.5 * ((p4[i] + p3[i]) - (p1[i] + p2[i]));
v[2][i] = p1[i] - p3[i];
}
return UnitCrossProd(v, nor);
}
inline int ProjectVector(double v[], const double n[])
{
// project 'v' on the plane with normal given by 'n' and then normalize 'v'
LinearCombination(InnerProd(n, n), v, -InnerProd(v, n), n, v);
return Normalize(v);
}
#endif // GLVIS_GEOM_UTILS_HPP
|
