collapse_degenerate.cpp 2.26 KB
Newer Older
Kirill Terekhov's avatar
Updates    
Kirill Terekhov 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
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
#include <stdio.h>
#include "inmost.h"


using namespace INMOST;
typedef Storage::real real;
typedef Storage::real_array real_array;

int main(int argc, char ** argv)
{
	if( argc < 2 )
	{
		printf("Usage: %s input_mesh [output_mesh]\n",argv[0]);
		return -1;
	}
	
	Mesh A("");
	A.Load(argv[1]);
	
	Mesh::GeomParam table;
	table[ORIENTATION] = FACE;
	table[BARYCENTER] = FACE | CELL | EDGE;
	A.PrepareGeometricData(table);
	
	rMatrix n(3,1), x(3,1), E(3,3), I(3,3), cx(3,1);
	I = rMatrix::Unit(3);
	
	double vol;
	int collapsed = 0;
	MarkerType rev = A.CreatePrivateMarker();
	for(Mesh::iteratorCell it = A.BeginCell(); it != A.EndCell(); ++it)
	{
		ElementArray<Face> faces = it->getFaces();
		// sum(x*n^T - (n^T*x)/3.0*I) = 0
		double vol = 0,vol0 = 0;
		E.Zero();
		it->Centroid(cx.data());
		A.FacesOrientation(faces,rev);
		
		for(ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); ++jt)
		{
			jt->Barycenter(x.data());
			jt->Normal(n.data());
			n*= jt->GetPrivateMarker(rev) ? -1.0:1.0;
			E += (x-cx)*n.Transpose();
			vol += n.DotProduct(x-cx)/3.0;
			vol0 += n.DotProduct(cx)/3.0;
		}
		if( vol < 0.0 )
		{
			vol = -vol;
			E = -E;
			vol0 = -vol0;
		}
		faces.RemPrivateMarker(rev);
		E = E-vol*I;
		if( E.FrobeniusNorm() > 1.0e-4 || vol0 > 1.0e-5)
		{
			std::cout << "Collapsing cell " << it->LocalID() << " err " << E.FrobeniusNorm() << std::endl;
			std::cout << "Volume " << it->Volume() << " vol " << vol << " vol0 " << vol0 << std::endl;
			std::cout << "E:" << std::endl;
			E.Print();
			//replace with node
			//it->Centroid(x.data());
			it->Delete();
			collapsed ++;
		}
	}
	A.ReleasePrivateMarker(rev);
	std::cout << "total collapsed: " << collapsed << std::endl;
	
	for(ElementType et = FACE; et > NONE; et = PrevElementType(et) )
		for(Mesh::iteratorElement it = A.BeginElement(et); it != A.EndElement(); ++it)
			if( it->nbAdjElements(CELL) == 0 ) it->Delete();
	
	
	if( A.HaveTag("GRIDNAME") )
	{
		Storage::bulk_array nameA = A.self().BulkArray(A.GetTag("GRIDNAME"));
		std::string ins = "_collapse_degenerate";
		nameA.insert(nameA.end(),ins.c_str(),ins.c_str()+ins.size());
	}
	
	if( argc > 2 )
	{
		std::cout << "Save to " << argv[2] << std::endl;
		A.Save(argv[2]);
	}
	else
	{
		std::cout << "Save to out.vtk" << std::endl;
		A.Save("out.vtk");
	}

	return 0;
}