Commit ad941614 authored by Heinrich Schindler's avatar Heinrich Schindler
Browse files

Reduce code duplication (functions GetElementVolume and GetAspectRatio).

parent 49df5686
Pipeline #3158095 passed with stage
in 6 minutes and 18 seconds
......@@ -44,12 +44,6 @@ class ComponentAnsys121 : public ComponentFieldMap {
/// Set the limits of the active region along z.
void SetRangeZ(const double zmin, const double zmax);
protected:
void UpdatePeriodicity() override;
double GetElementVolume(const unsigned int i) override;
void GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) override;
};
}
......
......@@ -36,13 +36,6 @@ class ComponentAnsys123 : public ComponentFieldMap {
bool SetWeightingField(std::string prnsol, std::string label);
protected:
// Verify periodicities
void UpdatePeriodicity() override { UpdatePeriodicityCommon(); }
double GetElementVolume(const unsigned int i) override;
void GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) override;
};
}
#endif
......@@ -197,11 +197,9 @@ class ComponentCST : public ComponentFieldMap {
unsigned int& i, unsigned int& j, unsigned int& k) const;
protected:
// Verify periodicities
void UpdatePeriodicity() override;
double GetElementVolume(const unsigned int i) override;
void GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) override;
double GetElementVolume(const size_t i) const override;
void GetAspectRatio(const size_t i,
double& dmin, double& dmax) const override;
/**
* Calculate the index in the vectors m_xlines, m_ylines, m_zlines, which is
......
......@@ -66,13 +66,6 @@ public:
void SetTimeInterval(const double mint, const double maxt,
const double stept);
protected:
void UpdatePeriodicity() override { UpdatePeriodicityCommon(); }
double GetElementVolume(const unsigned int i) override;
void GetAspectRatio(const unsigned int i, double &dmin,
double &dmax) override;
private:
double m_unit = 100.;
bool m_timeset = false;
......
// Copied and modified ComponentAnsys123.hh
#ifndef G_COMPONENT_ELMER_H
#define G_COMPONENT_ELMER_H
......@@ -52,13 +50,6 @@ class ComponentElmer : public ComponentFieldMap {
/// Import a list of voltages to be used as weighting field.
bool SetWeightingField(std::string prnsol, std::string label);
protected:
// Verify periodicities
void UpdatePeriodicity() override { UpdatePeriodicityCommon(); }
double GetElementVolume(const unsigned int i) override;
void GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) override;
};
}
#endif
......@@ -54,13 +54,6 @@ class ComponentElmer2D : public ComponentFieldMap {
void SetRangeZ(const double zmin, const double zmax);
protected:
// Verify periodicities
void UpdatePeriodicity() override;
double GetElementVolume(const unsigned int i) override;
void GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) override;
};
}
#endif
......@@ -65,7 +65,8 @@ class ComponentFieldMap : public Component {
/// Return the number of mesh elements.
virtual size_t GetNumberOfElements() const { return m_elements.size(); }
/// Return the volume and aspect ratio of a mesh element.
bool GetElement(const size_t i, double& vol, double& dmin, double& dmax);
bool GetElement(const size_t i, double& vol,
double& dmin, double& dmax) const;
/// Return the material and node indices of a mesh element.
virtual bool GetElement(const size_t i, size_t& mat, bool& drift,
std::vector<size_t>& nodes) const;
......@@ -187,6 +188,10 @@ class ComponentFieldMap : public Component {
void Prepare();
// Periodicities
void UpdatePeriodicity() override {
if (!m_is3d) UpdatePeriodicity2d();
UpdatePeriodicityCommon();
}
void UpdatePeriodicity2d();
void UpdatePeriodicityCommon();
......@@ -243,9 +248,9 @@ class ComponentFieldMap : public Component {
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,
double& dmax) = 0;
virtual double GetElementVolume(const size_t i) const;
virtual void GetAspectRatio(const size_t i,
double& dmin, double& dmax) const;
size_t GetWeightingFieldIndex(const std::string& label) const;
size_t GetOrCreateWeightingFieldIndex(const std::string& label);
......
......@@ -898,50 +898,4 @@ void ComponentAnsys121::SetRangeZ(const double zmin, const double zmax) {
m_maxBoundingBox[2] = m_mapmax[2] = std::max(zmin, zmax);
}
void ComponentAnsys121::UpdatePeriodicity() {
UpdatePeriodicity2d();
UpdatePeriodicityCommon();
}
double ComponentAnsys121::GetElementVolume(const unsigned int i) {
if (i >= m_elements.size()) return 0.;
const Element& element = m_elements[i];
const Node& n0 = m_nodes[element.emap[0]];
const Node& n1 = m_nodes[element.emap[1]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n3 = m_nodes[element.emap[3]];
const double surf =
0.5 *
(fabs((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)) +
fabs((n3.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n3.y - n0.y)));
return surf;
}
void ComponentAnsys121::GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) {
if (i >= m_elements.size()) {
dmin = dmax = 0.;
return;
}
const Element& element = m_elements[i];
const int np = 8;
// Loop over all pairs of vertices.
for (int j = 0; j < np - 1; ++j) {
const Node& nj = m_nodes[element.emap[j]];
for (int k = j + 1; k < np; ++k) {
const Node& nk = m_nodes[element.emap[k]];
// Compute distance.
const double dx = nj.x - nk.x;
const double dy = nj.y - nk.y;
const double dist = sqrt(dx * dx + dy * dy);
if (k == 1) {
dmin = dmax = dist;
} else {
if (dist < dmin) dmin = dist;
if (dist > dmax) dmax = dist;
}
}
}
}
}
......@@ -928,50 +928,4 @@ Medium* ComponentAnsys123::GetMedium(const double xin, const double yin,
return m_materials[element.matmap].medium;
}
double ComponentAnsys123::GetElementVolume(const unsigned int i) {
if (i >= m_elements.size()) return 0.;
const Element& element = m_elements[i];
const Node& n0 = m_nodes[element.emap[0]];
const Node& n1 = m_nodes[element.emap[1]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n3 = m_nodes[element.emap[3]];
const double vol =
fabs((n3.x - n0.x) *
((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
(n3.y - n0.y) *
((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
(n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
(n3.x - n0.x) * (n1.y - n0.y))) /
6.;
return vol;
}
void ComponentAnsys123::GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) {
if (i >= m_elements.size()) {
dmin = dmax = 0.;
return;
}
const Element& element = m_elements[i];
const int np = 4;
// Loop over all pairs of vertices.
for (int j = 0; j < np - 1; ++j) {
const Node& nj = m_nodes[element.emap[j]];
for (int k = j + 1; k < np; ++k) {
const Node& nk = m_nodes[element.emap[k]];
// Compute distance.
const double dx = nj.x - nk.x;
const double dy = nj.y - nk.y;
const double dz = nj.z - nk.z;
const double dist = sqrt(dx * dx + dy * dy + dz * dz);
if (k == 1) {
dmin = dmax = dist;
} else {
if (dist < dmin) dmin = dist;
if (dist > dmax) dmax = dist;
}
}
}
}
}
......@@ -1153,13 +1153,8 @@ int ComponentCST::Index2Element(const unsigned int i, const unsigned int j,
return i + j * (m_nx - 1) + k * (m_nx - 1) * (m_ny - 1);
}
void ComponentCST::UpdatePeriodicity() {
UpdatePeriodicity2d();
UpdatePeriodicityCommon();
}
void ComponentCST::GetAspectRatio(const unsigned int element, double& dmin,
double& dmax) {
void ComponentCST::GetAspectRatio(const size_t element, double& dmin,
double& dmax) const {
if (element >= m_nElements) {
dmin = dmax = 0.;
return;
......@@ -1173,7 +1168,7 @@ void ComponentCST::GetAspectRatio(const unsigned int element, double& dmin,
dmax = std::max({dx, dy, dz});
}
double ComponentCST::GetElementVolume(const unsigned int element) {
double ComponentCST::GetElementVolume(const size_t element) const {
if (element >= m_nElements) return 0.;
unsigned int i, j, k;
Element2Index(element, i, j, k);
......
......@@ -644,59 +644,6 @@ Medium *ComponentComsol::GetMedium(const double xin, const double yin,
return m_materials[element.matmap].medium;
}
double ComponentComsol::GetElementVolume(const unsigned int i) {
if (i >= m_elements.size())
return 0.;
const Element &element = m_elements[i];
const Node &n0 = m_nodes[element.emap[0]];
const Node &n1 = m_nodes[element.emap[1]];
const Node &n2 = m_nodes[element.emap[2]];
const Node &n3 = m_nodes[element.emap[3]];
// Uses formula V = |a (dot) b x c|/6
// with a => "3", b => "1", c => "2" and origin "0"
const double vol =
fabs((n3.x - n0.x) *
((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
(n3.y - n0.y) *
((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
(n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
(n3.x - n0.x) * (n1.y - n0.y))) /
6.;
return vol;
}
void ComponentComsol::GetAspectRatio(const unsigned int i, double &dmin,
double &dmax) {
if (i >= m_elements.size()) {
dmin = dmax = 0.;
return;
}
const Element &element = m_elements[i];
const int np = 4;
// Loop over all pairs of vertices.
for (int j = 0; j < np - 1; ++j) {
const Node &nj = m_nodes[element.emap[j]];
for (int k = j + 1; k < np; ++k) {
const Node &nk = m_nodes[element.emap[k]];
// Compute distance.
const double dx = nj.x - nk.x;
const double dy = nj.y - nk.y;
const double dz = nj.z - nk.z;
const double dist = sqrt(dx * dx + dy * dy + dz * dz);
if (k == 1) {
dmin = dmax = dist;
} else {
if (dist < dmin)
dmin = dist;
if (dist > dmax)
dmax = dist;
}
}
}
}
bool ComponentComsol::SetDelayedWeightingPotential(const std::string &field,
const std::string &label) {
if (!m_ready) {
......
......@@ -652,54 +652,4 @@ Medium* ComponentElmer::GetMedium(const double xin, const double yin,
return m_materials[element.matmap].medium;
}
double ComponentElmer::GetElementVolume(const unsigned int i) {
if (i >= m_elements.size()) return 0.;
const Element& element = m_elements[i];
const Node& n0 = m_nodes[element.emap[0]];
const Node& n1 = m_nodes[element.emap[1]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n3 = m_nodes[element.emap[3]];
// Uses formula V = |a (dot) b x c|/6
// with a => "3", b => "1", c => "2" and origin "0"
const double vol =
fabs((n3.x - n0.x) *
((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
(n3.y - n0.y) *
((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
(n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
(n3.x - n0.x) * (n1.y - n0.y))) /
6.;
return vol;
}
void ComponentElmer::GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) {
if (i >= m_elements.size()) {
dmin = dmax = 0.;
return;
}
const Element& element = m_elements[i];
const int np = 4;
// Loop over all pairs of vertices.
for (int j = 0; j < np - 1; ++j) {
const Node& nj = m_nodes[element.emap[j]];
for (int k = j + 1; k < np; ++k) {
const Node& nk = m_nodes[element.emap[k]];
// Compute distance.
const double dx = nj.x - nk.x;
const double dy = nj.y - nk.y;
const double dz = nj.z - nk.z;
const double dist = sqrt(dx * dx + dy * dy + dz * dz);
if (k == 1) {
dmin = dmax = dist;
} else {
if (dist < dmin) dmin = dist;
if (dist > dmax) dmax = dist;
}
}
}
}
} // namespace Garfield
......@@ -701,52 +701,4 @@ void ComponentElmer2D::SetRangeZ(const double zmin, const double zmax) {
m_maxBoundingBox[2] = m_mapmax[2] = std::max(zmin, zmax);
}
void ComponentElmer2D::UpdatePeriodicity() {
UpdatePeriodicity2d();
UpdatePeriodicityCommon();
}
double ComponentElmer2D::GetElementVolume(const unsigned int i) {
if (i >= m_elements.size()) return 0.;
const Element& element = m_elements[i];
const Node& n0 = m_nodes[element.emap[0]];
const Node& n1 = m_nodes[element.emap[1]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n3 = m_nodes[element.emap[3]];
const double surf =
0.5 *
(fabs((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)) +
fabs((n3.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n3.y - n0.y)));
return surf;
}
void ComponentElmer2D::GetAspectRatio(const unsigned int i, double& dmin,
double& dmax) {
if (i >= m_elements.size()) {
dmin = dmax = 0.;
return;
}
const Element& element = m_elements[i];
const int np = 8;
// Loop over all pairs of vertices.
for (int j = 0; j < np - 1; ++j) {
const Node& nj = m_nodes[element.emap[j]];
for (int k = j + 1; k < np; ++k) {
const Node& nk = m_nodes[element.emap[k]];
// Compute distance.
const double dx = nj.x - nk.x;
const double dy = nj.y - nk.y;
const double dist = sqrt(dx * dx + dy * dy);
if (k == 1) {
dmin = dmax = dist;
} else {
if (dist < dmin) dmin = dist;
if (dist > dmax) dmax = dist;
}
}
}
}
} // namespace Garfield
......@@ -122,7 +122,7 @@ Medium* ComponentFieldMap::GetMedium(const unsigned int imat) const {
}
bool ComponentFieldMap::GetElement(const size_t i, double& vol,
double& dmin, double& dmax) {
double& dmin, double& dmax) const {
if (i >= m_elements.size()) {
std::cerr << m_className << "::GetElement: Index out of range.\n";
return false;
......@@ -133,6 +133,91 @@ bool ComponentFieldMap::GetElement(const size_t i, double& vol,
return true;
}
double ComponentFieldMap::GetElementVolume(const size_t i) const {
if (i >= m_elements.size()) return 0.;
const Element& element = m_elements[i];
if (m_elementType == ElementType::CurvedTetrahedron) {
const Node& n0 = m_nodes[element.emap[0]];
const Node& n1 = m_nodes[element.emap[1]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n3 = m_nodes[element.emap[3]];
// Uses formula V = |a (dot) b x c|/6
// with a => "3", b => "1", c => "2" and origin "0"
const double vol =
fabs((n3.x - n0.x) * ((n1.y - n0.y) * (n2.z - n0.z) -
(n2.y - n0.y) * (n1.z - n0.z)) +
(n3.y - n0.y) * ((n1.z - n0.z) * (n2.x - n0.x) -
(n2.z - n0.z) * (n1.x - n0.x)) +
(n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
(n3.x - n0.x) * (n1.y - n0.y))) /
6.;
return vol;
} else if (m_elementType == ElementType::Serendipity) {
const Node& n0 = m_nodes[element.emap[0]];
const Node& n1 = m_nodes[element.emap[1]];
const Node& n2 = m_nodes[element.emap[2]];
const Node& n3 = m_nodes[element.emap[3]];
const double surf = 0.5 *
(fabs((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)) +
fabs((n3.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n3.y - n0.y)));
return surf;
}
return 0.;
}
void ComponentFieldMap::GetAspectRatio(const size_t i, double& dmin,
double& dmax) const {
if (i >= m_elements.size()) {
dmin = dmax = 0.;
return;
}
const Element& element = m_elements[i];
if (m_elementType == ElementType::CurvedTetrahedron) {
const int np = 4;
// Loop over all pairs of vertices.
for (int j = 0; j < np - 1; ++j) {
const Node& nj = m_nodes[element.emap[j]];
for (int k = j + 1; k < np; ++k) {
const Node& nk = m_nodes[element.emap[k]];
// Compute distance.
const double dx = nj.x - nk.x;
const double dy = nj.y - nk.y;
const double dz = nj.z - nk.z;
const double dist = sqrt(dx * dx + dy * dy + dz * dz);
if (k == 1) {
dmin = dmax = dist;
} else {
if (dist < dmin) dmin = dist;
if (dist > dmax) dmax = dist;
}
}
}
} else if (m_elementType == ElementType::Serendipity) {
const int np = 8;
// Loop over all pairs of vertices.
for (int j = 0; j < np - 1; ++j) {
const Node& nj = m_nodes[element.emap[j]];
for (int k = j + 1; k < np; ++k) {
const Node& nk = m_nodes[element.emap[k]];
// Compute distance.
const double dx = nj.x - nk.x;
const double dy = nj.y - nk.y;
const double dist = sqrt(dx * dx + dy * dy);
if (k == 1) {
dmin = dmax = dist;
} else {
if (dist < dmin) dmin = dist;
if (dist > dmax) dmax = dist;
}
}
}
}
}
bool ComponentFieldMap::GetElement(const size_t i, size_t& mat, bool& drift,
std::vector<size_t>& nodes) const {
if (i >= m_elements.size()) {
......
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