Skip to content

Commit

Permalink
Change Cartesian to AxisAligned, style fixes, update CHANGELOG
Browse files Browse the repository at this point in the history
  • Loading branch information
hughcars committed Aug 7, 2023
1 parent 8fca0b3 commit 6017be4
Show file tree
Hide file tree
Showing 6 changed files with 36 additions and 32 deletions.
4 changes: 3 additions & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,9 @@ The format of this changelog is based on
- Added Apptainer/Singularity container build definition for Palace.
- Added support for non axis aligned lumped ports and current sources. Key
words `"X"`, `"Y"`, `"Z"` and `"R"`, with optional prefix `"+"` or `"-"` still work,
but now alternative directions can be specified as vectors with 3 components.
but now alternative directions can be specified as vectors with 3
components. Users will be warned, an ultimately errored, if the specified
directions do not agree with axes directions discovered from the geometry.

## [0.11.2] - 2023-07-14

Expand Down
12 changes: 5 additions & 7 deletions examples/rings/mesh/mesh.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
using Gmsh: gmsh
using LinearAlgebra

function generate_ring_mesh(gui::Bool = false)
function generate_ring_mesh(gui::Bool=false)
kernel = gmsh.model.occ

gmsh.initialize()
Expand Down Expand Up @@ -113,8 +113,8 @@ function generate_ring_mesh(gui::Bool = false)
)

# Apply a rotation transformation to all entities in the model
rc = [0.0,0.0,0.0]
ra = [1.0,1.0,1.0]
rc = [0.0, 0.0, 0.0]
ra = [1.0, 1.0, 1.0]
θ = pi / 3
ra ./= norm(ra)
kernel.rotate(kernel.getEntities(), rc[1], rc[2], rc[3], ra[1], ra[2], ra[3], θ)
Expand All @@ -137,7 +137,6 @@ function generate_ring_mesh(gui::Bool = false)
inner_gap_group = gmsh.model.addPhysicalGroup(2, [inner_gap], -1, "hole_inner")
outer_gap_group = gmsh.model.addPhysicalGroup(2, [outer_gap], -1, "hole_outer")


# Generate mesh
gmsh.option.setNumber("Mesh.MeshSizeMin", l_ring)
gmsh.option.setNumber("Mesh.MeshSizeMax", l_farfield)
Expand Down Expand Up @@ -212,6 +211,5 @@ function generate_ring_mesh(gui::Bool = false)
gmsh.fltk.run()
end

gmsh.finalize()

end
return gmsh.finalize()
end
33 changes: 19 additions & 14 deletions palace/fem/lumpedelement.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@
#include <string>
#include <mfem.hpp>
#include "fem/integrator.hpp"
#include "utils/geodata.hpp"
#include "utils/communication.hpp"
#include "utils/geodata.hpp"

namespace palace
{
Expand Down Expand Up @@ -81,11 +81,8 @@ class UniformElementData : public LumpedElementData

// Check the user specified direction aligns with an axis direction
auto Dot = [](const auto &x_1, const auto &x_2)
{
return x_1[0] * x_2[0] + x_1[1] * x_2[1] + x_1[2] * x_2[2];
};
auto AbsAngle = [Dot](const auto &x_1, const auto &x_2)
{
{ return x_1[0] * x_2[0] + x_1[1] * x_2[1] + x_1[2] * x_2[2]; };
auto AbsAngle = [Dot](const auto &x_1, const auto &x_2) {
return std::acos(std::abs(Dot(x_1, x_2)) / std::sqrt(Dot(x_1, x_1) * Dot(x_2, x_2)));
};

Expand All @@ -97,21 +94,29 @@ class UniformElementData : public LumpedElementData
if (deviation_0_deg > angle_warning_deg && deviation_1_deg > angle_warning_deg)
{
auto normalized_0 = bounding_box.normals[0];
for (auto &x : normalized_0) { x /= std::sqrt(Dot(bounding_box.normals[0], bounding_box.normals[0])); }
for (auto &x : normalized_0)
{
x /= std::sqrt(Dot(bounding_box.normals[0], bounding_box.normals[0]));
}
auto normalized_1 = bounding_box.normals[1];
for (auto &x : normalized_1) { x /= std::sqrt(Dot(bounding_box.normals[1], bounding_box.normals[1])); }
Mpi::Warning("User specified direction {} does not align with either bounding box axis up to {} degrees.\n"
"Axis 1: {} ({} degrees)\nAxis 2: {} ({} degrees)",
input_dir, angle_warning_deg, normalized_0, deviation_0_deg, normalized_1, deviation_1_deg);
for (auto &x : normalized_1)
{
x /= std::sqrt(Dot(bounding_box.normals[1], bounding_box.normals[1]));
}
Mpi::Warning("User specified direction {} does not align with either bounding box "
"axis up to {} degrees.\n"
"Axis 1: {} ({} degrees)\nAxis 2: {} ({} degrees)",
input_dir, angle_warning_deg, normalized_0, deviation_0_deg,
normalized_1, deviation_1_deg);
}

MFEM_VERIFY(deviation_0_deg <= angle_error_deg || deviation_1_deg <= angle_error_deg,
"Specified direction does not align sufficiently with bounding box axes");
"Specified direction does not align sufficiently with bounding box axes");

// Compute the length from the most aligned normal direction.
l = 2 * std::sqrt(deviation_0_deg < deviation_1_deg
? Dot(bounding_box.normals[0], bounding_box.normals[0])
: Dot(bounding_box.normals[1], bounding_box.normals[1]));
? Dot(bounding_box.normals[0], bounding_box.normals[0])
: Dot(bounding_box.normals[1], bounding_box.normals[1]));
w = bounding_box.Area() / l;
}

Expand Down
1 change: 0 additions & 1 deletion palace/utils/configfile.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -815,7 +815,6 @@ auto ParseDataNode(json &j, const std::string &key_word, bool is_port = true)
}
}


return node;
}
} // namespace
Expand Down
10 changes: 5 additions & 5 deletions palace/utils/geodata.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -939,17 +939,17 @@ BoundingBall GetBoundingBall(mfem::ParMesh &mesh, int attr, bool bdr)
return GetBoundingBall(mesh, marker, bdr);
}

void GetCartesianBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
mfem::Vector &max)
void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
mfem::Vector &max)
{
mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
marker = 0;
marker[attr - 1] = 1;
GetCartesianBoundingBox(mesh, marker, bdr, min, max);
GetAxisAlignedBoundingBox(mesh, marker, bdr, min, max);
}

void GetCartesianBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
mfem::Vector &min, mfem::Vector &max)
void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
bool bdr, mfem::Vector &min, mfem::Vector &max)
{
int dim = mesh.SpaceDimension();
min.SetSize(dim);
Expand Down
8 changes: 4 additions & 4 deletions palace/utils/geodata.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,10 +52,10 @@ void AttrToMarker(int max_attr, const mfem::Array<int> &attrs, mfem::Array<int>
void AttrToMarker(int max_attr, const std::vector<int> &attrs, mfem::Array<int> &marker);

// Helper function to construct the bounding box for all elements with the given attribute.
void GetCartesianBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
mfem::Vector &max);
void GetCartesianBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
mfem::Vector &min, mfem::Vector &max);
void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
mfem::Vector &max);
void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
bool bdr, mfem::Vector &min, mfem::Vector &max);

// Struct describing a bounding box in terms of the center and face normals. The
// normals specify the distance from the center of the box.
Expand Down

0 comments on commit 6017be4

Please sign in to comment.