difference_same_match.cpp 9.06 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include "inmost.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

using namespace INMOST;


int main(int argc, char *argv[]) 
{
	if (argc < 3)
	{
		std::cout << "Usage: " << argv[0] ;
		std::cout << " mesh1 mesh2" << std::endl;
		return -1;
	}
	Mesh m1,m2;
	//m1.SetFileOption("VTK_GRID_DIMS","2");
	//m2.SetFileOption("VTK_GRID_DIMS","2");
20
21
	//m1.SetFileOption("VERBOSITY","2");
	//m2.SetFileOption("VERBOSITY","2");
22
23
24
25
26
27
28
29
30
	m1.Load(argv[1]);
	m2.Load(argv[2]);
	
	if( m1.NumberOfNodes() != m2.NumberOfNodes() )
	{
		std::cout << "Unequal number of nodes!" << std::endl;
		return -1;
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
31
	TagReal volumes_tag = m1.CreateTag("VOLUMES_FOR_NORMS_CALCULATION",DATA_REAL,CELL|FACE|EDGE|NODE,NONE,1);
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
	
	Tag ref12 = m1.CreateTag("M1TOM2",DATA_REMOTE_REFERENCE,NODE|EDGE|FACE|CELL,NONE,1);
	Tag ref21 = m2.CreateTag("M2TOM1",DATA_REMOTE_REFERENCE,NODE|EDGE|FACE|CELL,NONE,1);
	
	int nnodes = (int)m1.NumberOfNodes();
	
	
	
	std::vector<HandleType> nodes1(nnodes), nodes2(nnodes);
	int k;
	k = 0;
	for(Mesh::iteratorNode it = m1.BeginNode(); it != m1.EndNode(); ++it)
		nodes1[k++] = it->GetHandle();
	std::sort(nodes1.begin(),nodes1.end(),Mesh::CentroidComparator(&m1));
	k = 0;
	for(Mesh::iteratorNode it = m2.BeginNode(); it != m2.EndNode(); ++it)
		nodes2[k++] = it->GetHandle();
	std::sort(nodes2.begin(),nodes2.end(),Mesh::CentroidComparator(&m2));
	
	TagRealArray coords1(m1.CoordsTag()), coords2(m2.CoordsTag());
	
	int nmiss = 0;
	for(k = 0; k < nnodes; ++k)
	{
		if( (coords1(nodes1[k],3,1) - coords2(nodes2[k],3,1)).FrobeniusNorm() < 1.0e-8 )
		{
			m1.RemoteReference(nodes1[k],ref12) = RemoteHandleType(&m2,nodes2[k]);
			m2.RemoteReference(nodes2[k],ref21) = RemoteHandleType(&m1,nodes1[k]);
		}
		else
		{
			std::cout << "cannot match node " << k << " coords" << std::endl;
			std::cout << "at m1 "; coords1(nodes1[k],1,3).Print();
			std::cout << "at m2 "; coords1(nodes1[k],1,3).Print();
			nmiss++;
		}
	}
	if( nmiss )
		std::cout << "cannot match " << nmiss << " of NODE" << std::endl;
	std::cout << "matching done for NODE" << std::endl;
	
	std::vector<HandleType> remote_adj;
	
	for(ElementType etype = EDGE; etype <= CELL; etype = NextElementType(etype) )
	{
		nmiss = 0;
		for(Mesh::iteratorElement it = m1.BeginElement(etype); it != m1.EndElement(); ++it)
		{
			ElementArray<Element> adj = it->getAdjElements(PrevElementType(etype));
			for(ElementArray<Element>::iterator jt = adj.begin(); jt != adj.end(); ++jt)
				remote_adj.push_back(jt->RemoteReference(ref12).second);
			HandleType remote = m2.FindSharedAdjacency(&remote_adj[0],remote_adj.size());
			if( remote != InvalidHandle() )
				it->RemoteReference(ref12) = RemoteHandleType(&m2,remote);
			else
			{
				std::cout << "cannot match " << ElementTypeName(etype) << ":" << it->LocalID() << std::endl;
				nmiss++;
			}
			remote_adj.clear();
		}
		if( nmiss )
			std::cout << "cannot match " << nmiss << " of " << ElementTypeName(etype) << " from m1 to m2 " << std::endl;
		nmiss = 0;
		for(Mesh::iteratorElement it = m2.BeginElement(etype); it != m2.EndElement(); ++it)
		{
			ElementArray<Element> adj = it->getAdjElements(PrevElementType(etype));
			for(ElementArray<Element>::iterator jt = adj.begin(); jt != adj.end(); ++jt)
				remote_adj.push_back(jt->RemoteReference(ref21).second);
			HandleType remote = m1.FindSharedAdjacency(&remote_adj[0],remote_adj.size());
			if( remote != InvalidHandle() )
				it->RemoteReference(ref21) = RemoteHandleType(&m1,remote);
			else
			{
				std::cout << "cannot match " << ElementTypeName(etype) << ":" << it->LocalID() << std::endl;
				nmiss++;
			}
			remote_adj.clear();
		}
		if( nmiss )
			std::cout << "cannot match " << nmiss << " of " << ElementTypeName(etype) << " from m2 to m1 " << std::endl;
		std::cout << "matching done for " << ElementTypeName(etype) << std::endl;
	}
		
Kirill Terekhov's avatar
Kirill Terekhov committed
116
	/*
117
118
	for(ElementType etype = NODE; etype <= CELL; etype = etype << 1)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
119
120
121
#if defined(USE_OMP)
#pragma omp parallel for
#endif
122
123
124
125
126
127
128
129
130
131
132
133
134
135
		for(int id = 0; id < m1.LastLocalID(etype); id++)
		{
			Element c1 = m1.ElementByLocalID(etype,id);
			if( c1.isValid() )
			{
				Storage::real vol = 0;
				ElementArray<Cell> cells = c1.getCells();
				for(ElementArray<Cell>::iterator it = cells.begin(); it != cells.end(); ++it)
				{
					vol += it->Volume() / static_cast<Storage::real>(it->nbAdjElements(etype));
				}
				c1->RealDF(volumes_tag) = vol;
			}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
		std::cout << "volumes precomputed for " << ElementTypeName(etype) << std::endl;
	}
	*/
	for(int id = 0; id < m1.CellLastLocalID(); id++)
	{
		Cell c1 = m1.CellByLocalID(id);
		if( c1.isValid() )
		{
			double vol = c1.Volume();
			volumes_tag[c1] = vol;
			
			for(ElementType etype = NODE; etype <= FACE; etype = NextElementType(etype) )
			{
				ElementArray<Element> adj = c1.getAdjElements(etype);
				for(ElementArray<Element>::iterator it = adj.begin(); it != adj.end(); ++it)
					volumes_tag[*it] += vol / adj.size();
			}
		}
154
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
155
	std::cout << "volumes precomputed" << std::endl;
156
157
158
	
	for(Mesh::iteratorTag t = m1.BeginTag(); t != m1.EndTag(); t++)
	{
159
		//if( *t == m1.CoordsTag() ) continue;
Kirill Terekhov's avatar
Kirill Terekhov committed
160
		//if( t->GetSize() == ENUMUNDEF ) continue;
161
162
163
164
165
		if( t->GetDataType() != DATA_REAL 
#if defined(USE_AUTODIFF)
			&& t->GetDataType() != DATA_VARIABLE 
#endif
			) continue;
166
167
168
169
170
171
172
173
174
175
176
		if( m2.HaveTag(t->GetTagName()) )
		{
			Tag t2 = m2.GetTag(t->GetTagName());
			
			if( t->GetSize() != t2.GetSize() ) continue;
			
			for(ElementType etype = NODE; etype <= CELL; etype = NextElementType(etype))
			{
				if( m1.ElementDefined(*t,etype) && m2.ElementDefined(t2,etype) )
				{
					Storage::real Cnorm = 0, L1norm = 0, L2norm = 0, absval, Lvol = 0, vol;
177
					Storage::real max_val = -1.0e+20, min_val = 1.0e+20;
Kirill Terekhov's avatar
Kirill Terekhov committed
178
					int err = 0, tot = 0, diff_size = 0, tot_size = 0;
179
180
181
182
183
184
185
186
187
188
189
					for(int id = 0; id < m1.LastLocalID(etype); id++)
					{
						Element c1 = m1.ElementByLocalID(etype,id);
						Element c2 = MakeElement(c1.RemoteReference(ref12));
						if( c1.isValid() && c2.isValid() )
						{
							vol = c1->RealDF(volumes_tag);
							if( t->GetDataType() == DATA_REAL )
							{
								Storage::real_array arr1 = c1->RealArray(*t);
								Storage::real_array arr2 = c2->RealArray(t2);
Kirill Terekhov's avatar
Kirill Terekhov committed
190
191
192
193
194
								tot_size++;
								if( arr1.size() != arr2.size() )
								{
									diff_size++;
									continue;
195
									std::cout << "tag " << t->GetTagName() << " position " << k << "/" << arr1.size();
Kirill Terekhov's avatar
Kirill Terekhov committed
196
197
198
199
200
201
									std::cout << " arrays of different size ";
									std::cout << " on MESH1:" << ElementTypeName(etype) << ":" << c1.LocalID() << " " << arr1.size();
									std::cout << " on MESH2:" << ElementTypeName(etype) << ":" << c2.LocalID() << " " << arr2.size();
									std::cout << std::endl;
									continue;
								}
202
203
								for(int k = 0; k < (int)arr1.size(); k++)
								{
204
205
206
207
208
									max_val = std::max(max_val,arr1[k]);
									max_val = std::max(max_val,arr2[k]);
									min_val = std::min(min_val,arr1[k]);
									min_val = std::min(min_val,arr2[k]);

209
									absval = fabs(arr1[k]-arr2[k]);
Kirill Terekhov's avatar
Kirill Terekhov committed
210
										
211
									if( absval > 1.0e-3 )
Kirill Terekhov's avatar
Kirill Terekhov committed
212
									{
213
214
215
216
217
218
219
										if( true )
										{
											std::cout << "tag " << t->GetTagName() << " position " << k << "/" << arr1.size();
											std::cout << " on MESH1:" << ElementTypeName(etype) << ":" << c1.LocalID() << " " << arr1[k];
											std::cout << " on MESH2:" << ElementTypeName(etype) << ":" << c2.LocalID() << " " << arr2[k];
											std::cout << std::endl;
										}
Kirill Terekhov's avatar
Kirill Terekhov committed
220
										err++;
Kirill Terekhov's avatar
Kirill Terekhov committed
221
									}
222
									
Kirill Terekhov's avatar
Kirill Terekhov committed
223
224
225
226
									if( *t != m1.CoordsTag() )
										arr1[k] = absval;
									
									tot++;
227
228
229
230
231
232
233
									
									if( Cnorm < absval ) Cnorm = absval;
									L1norm += absval*vol;
									L2norm += absval*absval*vol;
									Lvol += vol;
								}
							}
234
#if defined(USE_AUTODIFF)
235
236
237
238
							else if( t->GetDataType() == DATA_VARIABLE )
							{
								Storage::var_array arr1 = c1->VariableArray(*t);
								Storage::var_array arr2 = c2->VariableArray(t2);
239
240
241
242
243
244
								tot_size++;
								if( arr1.size() != arr2.size() )
								{
									diff_size++;
									continue;
								}
245
246
								for(int k = 0; k < (int)arr1.size(); k++)
								{
247
248
249
250
251
252
									max_val = std::max(max_val,get_value(arr1[k]));
									max_val = std::max(max_val,get_value(arr2[k]));
									min_val = std::min(min_val,get_value(arr1[k]));
									min_val = std::min(min_val,get_value(arr2[k]));


253
254
255
256
									absval = fabs(get_value(arr1[k])-get_value(arr2[k]));
									
									arr1[k] = absval;
									
257
									if( absval > 1.0e-3 ) err++;
258
259
260
261
262
263
264
									
									if( Cnorm < absval ) Cnorm = absval;
									L1norm += absval*vol;
									L2norm += absval*absval*vol;
									Lvol += vol;
								}
							}
265
#endif
266
267
						}
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
268
269
					if( diff_size ) std::cout << "Size is different on " << diff_size << " / " << tot_size << " of " << ElementTypeName(etype) << std::endl;
					if( err ) std::cout << "Reported error on " << err << " / " << tot << " of " << ElementTypeName(etype) << std::endl;
270
					std::cout << "Data " << t->GetTagName() << " on " << ElementTypeName(etype) << " Cnorm " << Cnorm << " L1norm " << L1norm/Lvol << " L2norm " << sqrt(L2norm/Lvol) << " Min " << min_val << " Max " << max_val << std::endl;
271
272
273
274
				}
			}
		}
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
275
	m1.DeleteTag(ref12);
276
277
278
279
280
281
282
	m1.DeleteTag(volumes_tag);
	m1.Save("diff.gmv");
	m1.Save("diff.vtk");
	m1.Save("diff.pmf");
	
	return 0;
}