Skip to content

Facette scan #17

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
arowana755 wants to merge 1 commit into
base: master
Choose a base branch
from
+179 −72
Open

Facette scan #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
161 changes: 157 additions & 4 deletions src/mesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ double mesh::thiele(const int region /** region index, or -1 for all magnetic re
}


double mesh::avg(std::function<double(Nodes::Node, Nodes::index)> getter /**< [in] */,
double mesh::avg(const std::function<double(Nodes::Node, Nodes::index)>& getter /**< [in] */,
Nodes::index d /**< [in] */,
int region /**< region index, or -1 for all magnetic regions */) const
{
Expand All @@ -68,7 +68,7 @@ double mesh::avg(std::function<double(Nodes::Node, Nodes::index)> getter /**< [i
}

double mesh::doOnNodes(const double init_val, const Nodes::index coord,
std::function<bool(double, double)> whatToDo) const
const std::function<bool(double, double)>& whatToDo) const
{
double result(init_val);
std::for_each(node.begin(), node.end(),
Expand Down Expand Up @@ -100,7 +100,7 @@ void mesh::indexReorder()
std::for_each(fac.begin(), fac.end(), [this, &sf](Facette::Fac &fa)
{
int i0 = fa.ind[0], i1 = fa.ind[1], i2 = fa.ind[2];
std::set<Facette::Fac>::iterator it = sf.end();
auto it = sf.end();
for (int perm = 0; perm < 2; perm++)
{
for (int nrot = 0; nrot < 3; nrot++)
Expand Down Expand Up @@ -167,7 +167,107 @@ void mesh::sortNodes(Nodes::index long_axis)
});
}

double mesh::surface(std::vector<int> &facIndices)
void mesh::checkFacettes() const
{
std::vector<Facette::Fac> allFacCtnr;
Copy link
Copy Markdown
Contributor

@edgar-bonet edgar-bonet May 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The Facette::Fac class is too big for this usage: it has some vectors and matrices that are useful for mesh elements, but make no sense for the simple triangular faces of the tetrahedra. This is significant because we expect millions of such faces in the typical use case. A custom, leaner class/struc would be suitable here.

std::vector<std::pair<int, bool>> twoTetFacCtnr; // Bool: if the facettes are
// in the same region or not.
const int tetraN = 4;

// Put all facettes into allFacCtnr
std::for_each(tet.begin(), tet.end(), // For each tetrahedron
[this, &allFacCtnr, &twoTetFacCtnr](const Tetra::Tet &tetrahedron)
{
for (int i = 0; i < tetraN; i++) // For each 4 facettes
{
Facette::Fac curFac(node, 0, tetrahedron.idxPrm, {tetrahedron.ind[i],
tetrahedron.ind[(i + 1) % tetraN],
tetrahedron.ind[(i + 2) % tetraN]});

// Search for an equal facette, insert at the end if not found.
Copy link
Copy Markdown
Contributor

@edgar-bonet edgar-bonet May 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your comments are pretty nice: short, useful and to the point!

int j = 0;
while (j < allFacCtnr.size() && !(allFacCtnr[j] == curFac))
{ j++; }
Comment on lines +188 to +190
Copy link
Copy Markdown
Contributor

@edgar-bonet edgar-bonet May 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A for loop would be more idiomatic here.

More importantly, I am worried about the nesting of these loops. Looks to me like the algorithm is in O(n2), where n is the number of triangular faces. This can be annoying when n is in the millions. I do believe we can do way better.

if (j == allFacCtnr.size())
{ allFacCtnr.push_back(curFac); }
else
{ // If an equal facette is found
int k = 0;
while (k < twoTetFacCtnr.size() && twoTetFacCtnr[k].first != j)
{ k++; }

if (k < twoTetFacCtnr.size())
{ // If one facette belongs to 3 tetrahedrons, throws an error
std::cout << "Error : wrong mesh generation";
exit(1);
}
else
{ // If only 2 tetrahedrons, updates twoTetFacCtnr
twoTetFacCtnr.push_back({j, (curFac.idxPrm == allFacCtnr[j].idxPrm)});
}
}
}
});

std::sort(twoTetFacCtnr.begin(), twoTetFacCtnr.end());

int j = 0;
for (int i = 0; i < allFacCtnr.size(); i++)
{ // For each facette
if (j < twoTetFacCtnr.size() && i == twoTetFacCtnr[j].first && twoTetFacCtnr[j].second)
{ // If the current one belongs to two tetrahedrons of the same region
int k = 0;
while (k < fac.size() && !(fac[k] == allFacCtnr[i]))
{ k++; }
if (k < fac.size())
{
std::cout << "Error : an internal facette belonging "
"to a surface region has been found";
exit(1);
}
j++;
}
else if (j < twoTetFacCtnr.size() && i == twoTetFacCtnr[j].first && !twoTetFacCtnr[j].second)
{ // If the current one belongs to two tetrahedrons of differents regions
int k = 0;
while (k < fac.size() && !(fac[k] == allFacCtnr[i]))
{ k++; }
if (k < fac.size())
{ // Looking for a second identical facette
k++;
while (k < fac.size() && !(fac[k] == allFacCtnr[i]))
{ k++; }
if (k < fac.size())
{
std::cout << "Error : an interface facette belonging "
"to multiple surface regions has been found";
exit(1);
}
}
j++;
}
else
{ // If the current one belongs to only one tetrahedron
int k = 0;
while (k < fac.size() && !(fac[k] == allFacCtnr[i]))
{ k++; }
if (k < fac.size())
{ // Looking for a second identical facette
k++;
while (k < fac.size() && !(fac[k] == allFacCtnr[i]))
{ k++; }
if (k < fac.size())
{
std::cout << "Error : an surface facette belonging "
"to multiple surface regions has been found";
exit(1);
}
}
}
}
}

double mesh::surface(std::vector<int> &facIndices) const
{
double S(0);
std::for_each(facIndices.begin(),facIndices.end(),[this,&S](int idx)
Expand All @@ -190,3 +290,56 @@ void mesh::setExtSpaceField(Settings &s /**< [in] */)
k++;
});
}

void mesh::init_distrib(Settings const &mySets)
{
for (int nodeIdx = 0; nodeIdx < int(node.size()); ++nodeIdx)
{
Nodes::Node &n = node[nodeIdx];
n.d[Nodes::NEXT].phi = 0.;
n.d[Nodes::NEXT].phiv = 0.;
// A non-magnetic node's magnetization is a vector of NAN.
if (!magNode[nodeIdx])
{
n.d[Nodes::CURRENT].u = Eigen::Vector3d(NAN, NAN, NAN);
n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
continue;
}
// If the initial magnetization depends only on the node position, we do not need to
// build the list of region names.
if (mySets.getMagType() == POSITION_ONLY)
{
n.d[Nodes::CURRENT].u = mySets.getMagnetization(n.p);
n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
continue;
}
// Get the list of region indices this node belongs to.
std::set<int> nodeRegions;
// skip region 0 (__default__) which should be empty
for (size_t regIdx = 1; regIdx < volumeRegions.size(); ++regIdx)
{
const std::vector<int> &regionTetras = volumeRegions[regIdx];
bool node_in_region = std::any_of(EXEC_POL,
regionTetras.begin(), regionTetras.end(),
[this, nodeIdx](int tetIdx)
{
const Tetra::Tet &tetrahedron = tet[tetIdx];
for (int i = 0; i < Tetra::N; ++i) // node of tetrahedron tetIdx
{
if (tetrahedron.ind[i] == nodeIdx)
{ return true; }
}
return false;
});
if (node_in_region)
{ nodeRegions.insert(regIdx); }
}
// Get the list of region names this node belongs to.
std::vector<std::string> region_names;
region_names.resize(nodeRegions.size());
std::transform(nodeRegions.begin(), nodeRegions.end(), region_names.begin(),
[this](int regIdx){ return paramTetra[regIdx].regName; });
n.d[Nodes::CURRENT].u = mySets.getMagnetization(n.p, region_names);
n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
}
}
90 changes: 22 additions & 68 deletions src/mesh.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,8 @@ class mesh
{ magNode[te.ind[i]] = true; }
}
});

checkFacettes();

for(unsigned int i=0;i<fac.size();i++)
{
if (isMagnetic(fac[i]) && !mySets.paramFacette[fac[i].idxPrm].suppress_charges
Expand Down Expand Up @@ -177,7 +178,7 @@ class mesh
}

/** make_evol on i^th node */
inline void updateNode(int i, double vp, double vq, const double dt)
inline void updateNode(const int i, const double vp, const double vq, const double dt)
{ node[i].make_evol(vp*gamma0, vq*gamma0, dt); }

/** call evolution for all the nodes */
Expand All @@ -199,7 +200,7 @@ class mesh
/** total magnetic volume of the mesh */
double totalMagVol;

/** face container */
/** face container. Contains only those on a surface region. */
std::vector<Facette::Fac> fac;

/** tetrahedron container */
Expand All @@ -216,63 +217,7 @@ class mesh

/** computes an analytical initial magnetization distribution as a starting point for the
* simulation. If the node is not magnetic then it is set to NAN. */
inline void init_distrib(Settings const &mySets /**< [in] */)
{
for (int nodeIdx = 0; nodeIdx < int(node.size()); ++nodeIdx)
{
Nodes::Node &n = node[nodeIdx];
n.d[Nodes::NEXT].phi = 0.;
n.d[Nodes::NEXT].phiv = 0.;

// A non-magnetic node's magnetization is a vector of NAN.
if (!magNode[nodeIdx])
{
n.d[Nodes::CURRENT].u = Eigen::Vector3d(NAN, NAN, NAN);
n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
continue;
}

// If the initial magnetization depends only on the node position, we do not need to
// build the list of region names.
if (mySets.getMagType() == POSITION_ONLY)
{
n.d[Nodes::CURRENT].u = mySets.getMagnetization(n.p);
n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
continue;
}

// Get the list of region indices this node belongs to.
std::set<int> nodeRegions;
// skip region 0 (__default__) which should be empty
for (size_t regIdx = 1; regIdx < volumeRegions.size(); ++regIdx)
{
const std::vector<int> &regionTetras = volumeRegions[regIdx];
bool node_in_region = std::any_of(EXEC_POL,
regionTetras.begin(), regionTetras.end(),
[this, nodeIdx](int tetIdx)
{
const Tetra::Tet &tetrahedron = tet[tetIdx];
for (int i = 0; i < Tetra::N; ++i) // node of tetrahedron tetIdx
{
if (tetrahedron.ind[i] == nodeIdx)
{ return true; }
}
return false;
});
if (node_in_region)
{ nodeRegions.insert(regIdx); }
}

// Get the list of region names this node belongs to.
std::vector<std::string> region_names;
region_names.resize(nodeRegions.size());
std::transform(nodeRegions.begin(), nodeRegions.end(), region_names.begin(),
[this](int regIdx){ return paramTetra[regIdx].regName; });

n.d[Nodes::CURRENT].u = mySets.getMagnetization(n.p, region_names);
n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
}
}
void init_distrib(Settings const &mySets /**< [in] */);

/**
returns Thiele length in the magnetic volume region of index region.
Expand All @@ -285,7 +230,7 @@ class mesh
/**
average component of either u or v through getter on the whole set of tetetrahedron
*/
double avg(std::function<double(Nodes::Node, Nodes::index)> getter /**< [in] */,
double avg(const std::function<double(Nodes::Node, Nodes::index)>& getter /**< [in] */,
Nodes::index d /**< [in] */,
int region = -1 /**< region index, or -1 for all magnetic regions */) const;

Expand All @@ -294,7 +239,7 @@ class mesh
{
double min_dot_product = std::transform_reduce(EXEC_POL, edges.begin(), edges.end(), 1.0,
[](double a, double b){ return std::min(a, b); },
[this](const Edge edge)
[this](const Edge &edge)
{
Eigen::Vector3d m1 = getNode_u(edge.first);
Eigen::Vector3d m2 = getNode_u(edge.second);
Expand All @@ -318,15 +263,15 @@ class mesh
/** setter for node[i]; what_to_set will fix what is the part of the node struct to set (usefull
* for fmm_demag.h) */
inline void set(const int i /**< [in] */,
std::function<void(Nodes::Node &, const double)> what_to_set /**< [in] */,
const std::function<void(Nodes::Node &, const double)>& what_to_set /**< [in] */,
const double val /**< [in] */)
{ what_to_set(node[i], val); }

/** getter for the node_index value at position i */
inline int getNodeIndex(const int i) const { return node_index[i]; }

/** return true if this tetraedron is magnetic */
inline bool isMagnetic(const Tetra::Tet &t) { return (paramTetra[t.idxPrm].Ms > 0); }
inline bool isMagnetic(const Tetra::Tet &t) const { return (paramTetra[t.idxPrm].Ms > 0); }

/** return true if this facet is magnetic */
inline bool isMagnetic(const Facette::Fac &f)
Expand Down Expand Up @@ -387,7 +332,7 @@ class mesh
/** loop on nodes to apply predicate 'whatTodo' */
double doOnNodes(const double init_val /**< [in] */,
const Nodes::index coord /**< [in] */,
std::function<bool(double, double)> whatToDo /**< [in] */) const;
const std::function<bool(double, double)>& whatToDo /**< [in] */) const;

/** return the minimum of all nodes coordinate along coord axis */
inline double minNodes(const Nodes::index coord /**< [in] */) const
Expand All @@ -398,7 +343,8 @@ class mesh
/** return the maximum of all nodes coordinate along coord axis */
inline double maxNodes(const Nodes::index coord /**< [in] */) const
{
return doOnNodes(__DBL_MIN__, coord, [](double a, double b) { return a > b; });
return doOnNodes(__DBL_MIN__, coord, [](const double a, const double b)
{ return a > b; });
}

/** redefine orientation of triangular faces in accordance with the tetrahedron
Expand Down Expand Up @@ -430,14 +376,22 @@ class mesh
* the matrix we will have to solve for. */
void sortNodes(Nodes::index long_axis /**< [in] */);

/** Check if all facette verify one of the three following conditions :
* - Is part of 2 tetraedrons
* - Is part of 1 tetraedron and one surface region
* - Is part of 1 tetraedron and no surface region
* Throws an error if finds one which is not in one of those three cases.
*/
void checkFacettes() const;

/** returns the surface defined by the set of facets of indices in facIndices
* each elementary surface triangle defined by points p0,p1,p2 is computed using
* norm(cross(p0p1,p0p2))/2, it is always positive */
double surface(std::vector<int> &facIndices);
double surface(std::vector<int> &facIndices) const;

/**return true if facette f is already indexed in the list idxMagList.
* Uses operator== for Fac. */
bool isInMagList(std::vector<int> &idxMagList, Facette::Fac &f)
bool isInMagList(std::vector<int> &idxMagList, Facette::Fac &f) const
{
auto it = std::find_if(idxMagList.begin(),idxMagList.end(),
[this,&f](int idx) { return (fac[idx] == f); });
Expand Down