Commit 81ae57a9 authored by Heinrich Schindler's avatar Heinrich Schindler
Browse files

Merge branch 'fieldmaps' into 'master'

FEM field maps

See merge request !245
parents 6493b155 5ebaf538
Pipeline #3132174 passed with stage
in 5 minutes and 52 seconds
......@@ -75,7 +75,6 @@ class ComponentFieldMap : public Component {
// Options
void EnableCheckMapIndices(const bool on = true) {
m_checkMultipleElement = on;
if (on) m_lastElement = -1;
}
/// Option to eliminate mesh elements in conductors (default: on).
void EnableDeleteBackgroundElements(const bool on = true) {
......@@ -214,15 +213,15 @@ class ComponentFieldMap : public Component {
/// Find the element for a point in curved quadratic quadrilaterals.
int FindElement5(const double x, const double y, const double z, double& t1,
double& t2, double& t3, double& t4, double jac[4][4],
double& det);
double& det) const;
/// Find the element for a point in curved quadratic tetrahedra.
int FindElement13(const double x, const double y, const double z, double& t1,
double& t2, double& t3, double& t4, double jac[4][4],
double& det);
double& det) const;
/// Find the element for a point in a cube.
int FindElementCube(const double x, const double y, const double z,
double& t1, double& t2, double& t3, TMatrixD*& jac,
std::vector<TMatrixD*>& dN);
std::vector<TMatrixD*>& dN) const;
/// Move (xpos, ypos, zpos) to field map coordinates.
void MapCoordinates(double& xpos, double& ypos, double& zpos, bool& xmirrored,
......@@ -234,8 +233,8 @@ class ComponentFieldMap : public Component {
bool& zmirrored, double& rcoordinate,
double& rotation) const;
int ReadInteger(char* token, int def, bool& error);
double ReadDouble(char* token, double def, bool& error);
static int ReadInteger(char* token, int def, bool& error);
static double ReadDouble(char* token, double def, bool& error);
virtual double GetElementVolume(const unsigned int i) = 0;
virtual void GetAspectRatio(const unsigned int i, double& dmin,
......@@ -264,9 +263,6 @@ class ComponentFieldMap : public Component {
/// Flag to check if bounding boxes of elements are cached
bool m_cacheElemBoundingBoxes = false;
/// Keep track of the last element found.
int m_lastElement = -1;
/// Calculate local coordinates for curved quadratic triangles.
int Coordinates3(double x, double y, double z, double& t1, double& t2,
double& t3, double& t4, double jac[4][4], double& det,
......
......@@ -289,11 +289,7 @@ void ComponentFieldMap::Field13(const std::array<double, 10>& v,
int ComponentFieldMap::FindElement5(const double x, const double y,
double const z, double& t1, double& t2,
double& t3, double& t4, double jac[4][4],
double& det) {
// Check if bounding boxes of elements have been computed
if (!m_cacheElemBoundingBoxes) {
m_cacheElemBoundingBoxes = true;
}
double& det) const {
// Tetra list in the block that contains the input 3D point.
std::vector<int> tetList;
......@@ -308,22 +304,6 @@ int ComponentFieldMap::FindElement5(const double x, const double y,
// Initial values.
t1 = t2 = t3 = t4 = 0;
// Check previously used element
if (m_lastElement > -1 && !m_checkMultipleElement) {
const Element& element = m_elements[m_lastElement];
if (element.degenerate) {
if (Coordinates3(x, y, z, t1, t2, t3, t4, jac, det, element) == 0) {
if (t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 && t3 <= +1) {
return m_lastElement;
}
}
} else {
if (Coordinates5(x, y, z, t1, t2, t3, t4, jac, det, element) == 0) {
if (t1 >= -1 && t1 <= +1 && t2 >= -1 && t2 <= +1) return m_lastElement;
}
}
}
// Verify the count of volumes that contain the point.
int nfound = 0;
int imap = -1;
......@@ -347,7 +327,6 @@ int ComponentFieldMap::FindElement5(const double x, const double y,
if (t1 < 0 || t1 > 1 || t2 < 0 || t2 > 1 || t3 < 0 || t3 > 1) continue;
++nfound;
imap = idxToElemList;
m_lastElement = idxToElemList;
if (m_debug) {
std::cout << m_className << "::FindElement5:\n";
std::cout << " Found matching degenerate element " << idxToElemList
......@@ -374,7 +353,6 @@ int ComponentFieldMap::FindElement5(const double x, const double y,
if (t1 < -1 || t1 > 1 || t2 < -1 || t2 > 1) continue;
++nfound;
imap = idxToElemList;
m_lastElement = idxToElemList;
if (m_debug) {
std::cout << m_className << "::FindElement5:\n";
std::cout << " Found matching non-degenerate element "
......@@ -404,7 +382,6 @@ int ComponentFieldMap::FindElement5(const double x, const double y,
std::cout << " No element matching point (" << x << ", " << y
<< ") found.\n";
}
m_lastElement = -1;
return -1;
}
if (nfound > 1) {
......@@ -421,7 +398,6 @@ int ComponentFieldMap::FindElement5(const double x, const double y,
t3 = t3bak;
t4 = t4bak;
imap = imapbak;
m_lastElement = imap;
return imap;
}
......@@ -436,7 +412,7 @@ int ComponentFieldMap::FindElement5(const double x, const double y,
int ComponentFieldMap::FindElement13(const double x, const double y,
const double z, double& t1, double& t2,
double& t3, double& t4, double jac[4][4],
double& det) {
double& det) const {
// Backup
double jacbak[4][4];
double detbak = 1.;
......@@ -446,17 +422,6 @@ int ComponentFieldMap::FindElement13(const double x, const double y,
// Initial values.
t1 = t2 = t3 = t4 = 0.;
// Check previously used element
if (m_lastElement > -1 && !m_checkMultipleElement) {
const Element& element = m_elements[m_lastElement];
if (Coordinates13(x, y, z, t1, t2, t3, t4, jac, det, element) == 0) {
if (t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 && t3 <= +1 &&
t4 >= 0 && t4 <= +1) {
return m_lastElement;
}
}
}
// Tetra list in the block that contains the input 3D point.
std::vector<int> tetList;
if (m_useTetrahedralTree && m_octree) {
......@@ -487,7 +452,6 @@ int ComponentFieldMap::FindElement13(const double x, const double y,
}
++nfound;
imap = idxToElemList;
m_lastElement = idxToElemList;
if (m_debug) {
std::cout << m_className << "::FindElement13:\n";
std::cout << " Found matching element " << i << ".\n";
......@@ -515,7 +479,6 @@ int ComponentFieldMap::FindElement13(const double x, const double y,
std::cout << " No element matching point (" << x << ", " << y << ", "
<< z << ") found.\n";
}
m_lastElement = -1;
return -1;
}
if (nfound > 1) {
......@@ -532,7 +495,6 @@ int ComponentFieldMap::FindElement13(const double x, const double y,
t3 = t3bak;
t4 = t4bak;
imap = imapbak;
m_lastElement = imap;
return imap;
}
......@@ -547,42 +509,26 @@ int ComponentFieldMap::FindElement13(const double x, const double y,
int ComponentFieldMap::FindElementCube(const double x, const double y,
const double z, double& t1, double& t2,
double& t3, TMatrixD*& jac,
std::vector<TMatrixD*>& dN) {
std::vector<TMatrixD*>& dN) const {
int imap = -1;
if (m_lastElement >= 0) {
const Element& element = m_elements[m_lastElement];
const size_t nElements = m_elements.size();
for (size_t i = 0; i < nElements; ++i) {
const Element& element = m_elements[i];
const Node& n3 = m_nodes[element.emap[3]];
if (x >= n3.x && y >= n3.y && z >= n3.z) {
const Node& n0 = m_nodes[element.emap[0]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n7 = m_nodes[element.emap[7]];
if (x < n0.x && y < n2.y && z < n7.z) {
imap = m_lastElement;
}
}
}
// Default element loop
if (imap == -1) {
const size_t nElements = m_elements.size();
for (size_t i = 0; i < nElements; ++i) {
const Element& element = m_elements[i];
const Node& n3 = m_nodes[element.emap[3]];
if (x < n3.x || y < n3.y || z < n3.z) continue;
const Node& n0 = m_nodes[element.emap[0]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n7 = m_nodes[element.emap[7]];
if (x < n0.x && y < n2.y && z < n7.z) {
imap = i;
break;
}
if (x < n3.x || y < n3.y || z < n3.z) continue;
const Node& n0 = m_nodes[element.emap[0]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n7 = m_nodes[element.emap[7]];
if (x < n0.x && y < n2.y && z < n7.z) {
imap = i;
break;
}
}
if (imap < 0) {
if (m_debug) {
std::cout << m_className << "::FindElementCube:\n";
std::cout << " Point (" << x << "," << y << "," << z
std::cout << m_className << "::FindElementCube:\n"
<< " Point (" << x << "," << y << "," << z
<< ") not in the mesh, it is background or PEC.\n";
const Node& first0 = m_nodes[m_elements.front().emap[0]];
const Node& first2 = m_nodes[m_elements.front().emap[2]];
......@@ -862,8 +808,8 @@ int ComponentFieldMap::Coordinates3(const double x, const double y,
double& t3, double& t4, double jac[4][4],
double& det, const Element& element) const {
if (m_debug) {
std::cout << m_className << "::Coordinates3:\n";
std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
std::cout << m_className << "::Coordinates3:\n"
<< " Point (" << x << ", " << y << ", " << z << ")\n";
}
// Provisional values
......@@ -880,8 +826,8 @@ int ComponentFieldMap::Coordinates3(const double x, const double y,
const double d3 =
(n2.x - n0.x) * (n1.y - n0.y) - (n1.x - n0.x) * (n2.y - n0.y);
if (d1 == 0 || d2 == 0 || d3 == 0) {
std::cerr << m_className << "::Coordinates3:\n";
std::cerr << " Calculation of linear coordinates failed; abandoned.\n";
std::cerr << m_className << "::Coordinates3:\n"
<< " Calculation of linear coordinates failed; abandoned.\n";
return 1;
}
t1 = ((x - n1.x) * (n2.y - n1.y) - (y - n1.y) * (n2.x - n1.x)) / d1;
......@@ -1007,8 +953,8 @@ int ComponentFieldMap::Coordinates4(const double x, const double y,
const Element& element) const {
// Debugging
if (m_debug) {
std::cout << m_className << "::Coordinates4:\n";
std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
std::cout << m_className << "::Coordinates4:\n"
<< " Point (" << x << ", " << y << ", " << z << ")\n";
}
// Failure flag
......@@ -1382,8 +1328,8 @@ int ComponentFieldMap::Coordinates13(const double x, const double y,
double& det,
const Element& element) const {
if (m_debug) {
std::cout << m_className << "::Coordinates13:\n";
std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
std::cout << m_className << "::Coordinates13:\n"
<< " Point (" << x << ", " << y << ", " << z << ")\n";
}
// Provisional values
......@@ -1397,9 +1343,9 @@ int ComponentFieldMap::Coordinates13(const double x, const double y,
if (t1 < -f || t2 < -f || t3 < -f || t4 < -f || t1 > 1 + f || t2 > 1 + f ||
t3 > 1 + f || t4 > 1 + f) {
if (m_debug) {
std::cout << m_className << "::Coordinates13:\n";
std::cout << " Linear isoparametric coordinates more than\n";
std::cout << " f (" << f << ") out of range.\n";
std::cout << m_className << "::Coordinates13:\n"
<< " Linear isoparametric coordinates more than\n"
<< " f (" << f << ") out of range.\n";
}
return 0;
}
......@@ -1646,7 +1592,6 @@ void ComponentFieldMap::Reset() {
m_octree.reset(nullptr);
m_cacheElemBoundingBoxes = false;
m_lastElement = -1;
}
void ComponentFieldMap::Prepare() {
......
......@@ -73,6 +73,11 @@ bool ViewFEMesh::Plot(const bool twod) {
pad->cd();
if (!twod) {
if (!m_cmp->m_is3d) {
std::cerr << m_className << "::Plot:\n"
<< " Cannot plot 2D mesh elements in 3D.\n";
return false;
}
DrawElements3d();
DrawDriftLines3d();
return true;
......
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