HexMeshLib/JacobianEva.cpp at master · yuwuyi/HexMeshLib · GitHub

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
#include "JacobianEva.h"
#include "HexVertex.h"
#include "Hex.h"
#include "HexMesh.h"
#include "HexIterators.h"
#include "Eigen/Eigen"
#include <numeric>
#include <iostream>

using namespace Eigen;

JacobianEva::JacobianEva(HexMesh *hm) : m_hexmesh (hm) {

}

JacobianEva::~JacobianEva() {

}

double JacobianEva::evaluate(HexVertex *hv, HexVertex* nhv[3]) {
	Matrix3d J = Matrix3d::Zero();
	for (int i = 0; i < 3; ++i) {
		Point pt = nhv[i]->point() - hv->point();
		pt /= pt.norm();
		for (int j = 0; j < 3; ++j) {
			J(j, i) = pt[j];
		}
	}

	return J.determinant();
}

void JacobianEva::evaluate(Hex *hex) {
	int corr[8][3] = {
		{4,2,1},
		{0,3,5},
		{0,6,3},
		{1,2,7},
		{6,0,5},
		{1,7,4},
		{2,4,7},
		{6,5,3},
	};

	double min_q = 1e10;

	for (int i = 0; i < 8; ++i) {
		HexVertex *hvs[3] = {NULL, NULL, NULL};
		for (int j = 0; j < 3; ++j) {
			hvs[j] = hex->vertex(corr[i][j]);
		}
		double q = evaluate(hex->vertex(i), hvs);

		if (q < min_q) {
			min_q = q;
		}
		vertexmap[hex->vertex(i)].push_back(q);
	}

	return;
}

void JacobianEva::evaluate() {
	vertexmap.clear();

	for (HexMeshHexIterator hmhit(m_hexmesh); !hmhit.end(); ++hmhit)  {
		Hex *hex = *hmhit;
		evaluate(hex);
	}

	for (std::map<HexVertex*, std::vector<double> >::iterator b = vertexmap.begin(), e = vertexmap.end(); b != e; ++b) {
		HexVertex *hv = (*b).first;
		std::vector<double> &vec = (*b).second;
		vertexJacobian[hv] = std::accumulate(vec.begin(), vec.end(), 0.0) / vec.size();
	}

}