Commit 0096fa49 authored by Kirill Terekhov's avatar Kirill Terekhov

Mesh data difference with different order of elements

Allows to compare results on meshes obtained from parallel runs with
different number of processors
parent 2dce79a5
Pipeline #222 failed with stages
in 10 minutes and 23 seconds
......@@ -16,7 +16,8 @@ add_executable(Sew sew.cpp)
add_executable(Scale scale.cpp)
add_executable(Move move.cpp)
add_executable(Convert convert.cpp)
add_executable(SameMeshDifference difference_same.cpp)
add_executable(SameMeshDifference difference_same.cpp)
add_executable(MatchSameMeshDifference difference_same_match.cpp)
add_executable(MeshDifference difference_map.cpp)
add_executable(Kmeans kmeans.cpp)
add_executable(Agglomerate agglomerate.cpp)
......@@ -268,7 +269,18 @@ if(USE_MPI)
set_target_properties(SameMeshDifference PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS SameMeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS SameMeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(MatchSameMeshDifference inmost)
if(USE_MPI)
message("linking MatchSameMeshDifference with MPI")
target_link_libraries(MatchSameMeshDifference ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(MatchSameMeshDifference PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS MatchSameMeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(MeshDifference inmost)
......
#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");
m1.SetFileOption("VERBOSITY","2");
m2.SetFileOption("VERBOSITY","2");
m1.Load(argv[1]);
m2.Load(argv[2]);
if( m1.NumberOfNodes() != m2.NumberOfNodes() )
{
std::cout << "Unequal number of nodes!" << std::endl;
return -1;
}
Tag volumes_tag = m1.CreateTag("VOLUMES_FOR_NORMS_CALCULATION",DATA_REAL,CELL|FACE|EDGE|NODE,NONE,1);
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;
}
for(ElementType etype = NODE; etype <= CELL; etype = etype << 1)
{
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;
}
}
}
for(Mesh::iteratorTag t = m1.BeginTag(); t != m1.EndTag(); t++)
{
if( *t == m1.CoordsTag() ) continue;
if( t->GetSize() == ENUMUNDEF ) continue;
if( t->GetDataType() != DATA_REAL && t->GetDataType() != DATA_VARIABLE ) continue;
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;
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);
for(int k = 0; k < (int)arr1.size(); k++)
{
absval = fabs(arr1[k]-arr2[k]);
arr1[k] = absval;
if( Cnorm < absval ) Cnorm = absval;
L1norm += absval*vol;
L2norm += absval*absval*vol;
Lvol += vol;
}
}
else if( t->GetDataType() == DATA_VARIABLE )
{
Storage::var_array arr1 = c1->VariableArray(*t);
Storage::var_array arr2 = c2->VariableArray(t2);
for(int k = 0; k < (int)arr1.size(); k++)
{
absval = fabs(get_value(arr1[k])-get_value(arr2[k]));
arr1[k] = absval;
if( Cnorm < absval ) Cnorm = absval;
L1norm += absval*vol;
L2norm += absval*absval*vol;
Lvol += vol;
}
}
}
}
std::cout << "Data " << t->GetTagName() << " on " << ElementTypeName(etype) << " Cnorm " << Cnorm << " L1norm " << L1norm/Lvol << " L2norm " << sqrt(L2norm/Lvol) << std::endl;
}
}
}
}
m1.DeleteTag(volumes_tag);
m1.Save("diff.gmv");
m1.Save("diff.vtk");
m1.Save("diff.pmf");
return 0;
}
......@@ -32,6 +32,20 @@ difference_same mesh_file1 mesh_file2
mesh_file1 - first mesh file
mesh_file2 - second mesh file
difference_same_match - Computes difference between the real data on all elements of two similar grids
and writes the grid with the difference. It writes absolute value of difference.
These variant works for meshes that has different order of elements but same mesh
geometry.
Outputs L_inf, L_1, L_2 norms of the difference.
Similar grids means their geometry should match between each other
Usage:
difference_same_match mesh_file1 mesh_file2
mesh_file1 - first mesh file
mesh_file2 - second mesh file
difference_map - Computes difference between the real data on cells of two different grids
and writes both grids with the difference. It writes absolute value of difference.
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment