Write mesh to output file.
{
int surf1 =
m_config[
"surf1"]. as<int> ();
int surf2 =
m_config[
"surf2"]. as<int> ();
string dir =
m_config[
"dir"]. as<string>();
bool orient =
m_config[
"orient"].as<
bool> ();
if (surf1 == -1)
{
cerr << "WARNING: surf1 must be set to a positive integer. "
<< "Skipping periodic alignment." << endl;
return;
}
if (surf2 == -1)
{
cerr << "WARNING: surf2 must be set to a positive integer. "
<< "Skipping periodic alignment." << endl;
return;
}
if (dir != "x" && dir != "y" && dir != "z")
{
cerr << "WARNING: dir must be set to either x, y or z. "
<< "Skipping periodic alignment." << endl;
return;
}
vec[0] = dir == "x" ? 1.0 : 0.0;
vec[1] = dir == "y" ? 1.0 : 0.0;
vec[2] = dir == "z" ? 1.0 : 0.0;
if (it1 ==
m_mesh->m_composite.end())
{
cerr << "WARNING: Couldn't find surface " << surf1
<< ". Skipping periodic alignment." << endl;
return;
}
if (it2 ==
m_mesh->m_composite.end())
{
cerr << "WARNING: Couldn't find surface " << surf2 << ", "
<< "skipping periodic alignment." << endl;
return;
}
if (c1->m_items.size() != c2->m_items.size())
{
cerr << "WARNING: Surfaces " << surf1 << " and " << surf2
<< " have different numbers of elements. Skipping periodic"
<< " alignment." << endl;
return;
}
c1->m_reorder = false;
c2->m_reorder = false;
map<int, pair<FaceSharedPtr, vector<int> > > perFaces;
map<int, Node> centroidMap;
for (int i = 0; i < c2->m_items.size(); ++i)
{
Node centroid;
for (int j = 0; j < c2->m_items[i]->GetVertexCount(); ++j)
{
centroid += *(c2->m_items[i]->GetVertex(j));
}
centroid /= (
NekDouble)c2->m_items[i]->GetVertexCount();
centroidMap[i] = centroid;
}
boost::unordered_set<int> elmtDone;
map<int, int> elmtPairs;
map<int, int> vertCheck;
for (int i = 0; i < c1->m_items.size(); ++i)
{
Node centroid;
for (int j = 0; j < c1->m_items[i]->GetVertexCount(); ++j)
{
centroid += *(c1->m_items[i]->GetVertex(j));
}
centroid /= (
NekDouble)c1->m_items[i]->GetVertexCount();
for (it = centroidMap.begin(); it != centroidMap.end(); ++it)
{
if (elmtDone.count(it->first) > 0)
{
continue;
}
Node dx = it->second - centroid;
if (fabs(fabs(dx.m_x*vec[0] + dx.m_y*vec[1] + dx.m_z*vec[2])/
sqrt(dx.abs2()) - 1.0) < 1e-8)
{
int id1, id2;
{
id1 = c1->m_items[i] ->GetEdgeLink()->m_id;
id2 = c2->m_items[it->first]->GetEdgeLink()->m_id;
}
else
{
id1 = c1->m_items[i] ->GetFaceLink()->m_id;
id2 = c2->m_items[it->first]->GetFaceLink()->m_id;
}
elmtDone.insert(it->first);
elmtPairs[i] = it->first;
int nVerts = c1->m_items[i]->GetVertexCount();
vector<int> perVerts(nVerts, 0), perVertsInv(nVerts, 0);
if (orient)
{
for (int k = 0; k < nVerts; ++k)
{
NodeSharedPtr n1 = c1->m_items[i]->GetFaceLink()->m_vertexList[k];
int l;
for (l = 0; l < nVerts; ++l)
{
c2->m_items[it->first]->GetFaceLink()->m_vertexList[l];
Node dn = *n2 - *n1;
if (fabs(fabs(dn.m_x*vec[0] + dn.m_y*vec[1] +
dn.m_z*vec[2])/
sqrt(dn.abs2()) - 1.0) < 1e-8)
{
perVerts [k] = l;
perVertsInv[l] = k;
int id1 = n1->m_id;
int id2 = n2->m_id;
if (vertCheck.count(id1) == 0)
{
vertCheck[id1] = id2;
}
else
{
"Periodic vertex already "
"identified!");
}
break;
}
}
"Could not identify periodic vertices.");
}
int tot1 = 0, tot2 = 0;
for (int k = 0; k < nVerts; ++k)
{
tot1 += perVerts [k];
tot2 += perVertsInv[k];
}
tot2 == nVerts*(nVerts-1)/2,
"Error identifying periodic vertices");
}
{
perFaces[id1] = make_pair(
c2->m_items[it->first]->GetFaceLink(), perVerts);
perFaces[id2] = make_pair(
c1->m_items[i] ->GetFaceLink(), perVertsInv);
}
break;
}
}
if (it == centroidMap.end())
{
cerr << "WARNING: Could not find matching edge for surface "
<< "element " << c1->m_items[i]->GetId() << ". "
<< "Skipping periodic alignment." << endl;
return;
}
}
vector<ElementSharedPtr> tmp = c2->m_items;
for (int i = 0; i < tmp.size(); ++i)
{
c2->m_items[i] = tmp[elmtPairs[i]];
}
if (orient)
{
}
}