diff --git a/.gitignore b/.gitignore index 49ab51382..9cf12d78f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ build*/ examples/**/postpro*/ examples/**/log*/ +examples/**/*.log spack/local/packages/palace/__pycache__/ test/ref/**/*.json .*.swp diff --git a/CHANGELOG.md b/CHANGELOG.md index b1613f445..2a784f99b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,12 @@ The format of this changelog is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/). +## [In progress] + + - 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. + ## [0.11.2] - 2023-07-14 - Changed layout and names of `palace/` source directory for better organization. diff --git a/examples/cpw/cpw_lumped_uniform.json b/examples/cpw/cpw_lumped_uniform.json index e3bce2129..e5e62cd43 100644 --- a/examples/cpw/cpw_lumped_uniform.json +++ b/examples/cpw/cpw_lumped_uniform.json @@ -79,11 +79,11 @@ [ { "Attributes": [4], - "Direction": "+Y" + "Direction": [0,1,0] }, { "Attributes": [8], - "Direction": "-Y" + "Direction": [0,-1,0] } ] }, @@ -94,11 +94,11 @@ [ { "Attributes": [5], - "Direction": "+Y" + "Direction": [0,1,0] }, { "Attributes": [9], - "Direction": "-Y" + "Direction": [0,-1,0] } ] }, @@ -109,11 +109,11 @@ [ { "Attributes": [6], - "Direction": "+Y" + "Direction": [0,1,0] }, { "Attributes": [10], - "Direction": "-Y" + "Direction": [0,-1,0] } ] }, @@ -124,11 +124,11 @@ [ { "Attributes": [7], - "Direction": "+Y" + "Direction": [0,1,0] }, { "Attributes": [11], - "Direction": "-Y" + "Direction": [0,-1,0] } ] } diff --git a/palace/fem/lumpedelement.hpp b/palace/fem/lumpedelement.hpp index 4e98a6e81..50a9b850e 100644 --- a/palace/fem/lumpedelement.hpp +++ b/palace/fem/lumpedelement.hpp @@ -55,38 +55,22 @@ class LumpedElementData class UniformElementData : public LumpedElementData { protected: - bool sign; // Sign of incident field, +x̂ / ŷ / ẑ if true - int component; // Lumped element direction (0: x, 1: y, 2: z) - double l, w; // Lumped element length and width + mfem::Vector direction; // Cartesian vector specifying signed direction of incident field + double l, w; // Lumped element length and width public: - UniformElementData(const std::string &direction, const mfem::Array &marker, + UniformElementData(const std::array &input_dir, const mfem::Array &marker, mfem::ParFiniteElementSpace &fespace) - : LumpedElementData(fespace.GetParMesh()->SpaceDimension(), marker), - sign(direction[0] == '+') + : LumpedElementData(fespace.GetParMesh()->SpaceDimension(), marker), direction(3) { - switch (direction[1]) - { - case 'x': - component = 0; - break; - case 'y': - component = 1; - break; - case 'z': - component = 2; - break; - default: - MFEM_ABORT("Lumped element direction is not correctly formatted!"); - component = 0; // For compiler warning - break; - } + std::copy(input_dir.begin(), input_dir.end(), direction.begin()); - // Get the lumped element length and width assuming axis-aligned rectangle. + // Get the lumped element length and width. mfem::Vector bbmin, bbmax; mesh::GetBoundingBox(*fespace.GetParMesh(), marker, true, bbmin, bbmax); double A = GetArea(fespace); - l = bbmax(component) - bbmin(component); + bbmax -= bbmin; + l = std::abs(bbmax * direction); w = A / l; } @@ -96,9 +80,8 @@ class UniformElementData : public LumpedElementData std::unique_ptr GetModeCoefficient(double coef = 1.0) const override { - mfem::Vector source(dim); - source = 0.0; - source(component) = (sign ? 1.0 : -1.0) * coef; + mfem::Vector source = direction; + source *= coef; return std::make_unique(source); } }; @@ -151,7 +134,7 @@ class CoaxialElementData : public LumpedElementData { double scoef = (sign ? 1.0 : -1.0) * coef; mfem::Vector x0(c); - auto Source = [scoef, x0](const mfem::Vector &x, mfem::Vector &f) -> void + auto Source = [scoef, x0](const mfem::Vector &x, mfem::Vector &f) { f = x; f -= x0; diff --git a/palace/models/lumpedportoperator.cpp b/palace/models/lumpedportoperator.cpp index acbd9717d..8d752a487 100644 --- a/palace/models/lumpedportoperator.cpp +++ b/palace/models/lumpedportoperator.cpp @@ -50,13 +50,6 @@ LumpedPortData::LumpedPortData(const config::LumpedPortData &data, // Construct the port elements allowing for a possible multielement lumped port. for (const auto &node : data.nodes) { - // Check input direction. - MFEM_VERIFY(node.direction.length() == 2 && - (node.direction[0] == '-' || node.direction[0] == '+') && - (node.direction[1] == 'x' || node.direction[1] == 'y' || - node.direction[1] == 'z' || node.direction[1] == 'r'), - "Lumped port direction is not correctly formatted!"); - mfem::Array attr_marker; mesh::AttrToMarker(h1_fespace.GetParMesh()->bdr_attributes.Max(), node.attributes, attr_marker); @@ -68,7 +61,7 @@ LumpedPortData::LumpedPortData(const config::LumpedPortData &data, else { elems.push_back( - std::make_unique(node.direction, attr_marker, h1_fespace)); + std::make_unique(node.normal, attr_marker, h1_fespace)); } } diff --git a/palace/models/surfacecurrentoperator.cpp b/palace/models/surfacecurrentoperator.cpp index 507b889d1..2a0d8b4f1 100644 --- a/palace/models/surfacecurrentoperator.cpp +++ b/palace/models/surfacecurrentoperator.cpp @@ -19,13 +19,6 @@ SurfaceCurrentData::SurfaceCurrentData(const config::SurfaceCurrentData &data, // sources. for (const auto &node : data.nodes) { - // Check input direction. - MFEM_VERIFY(node.direction.length() == 2 && - (node.direction[0] == '-' || node.direction[0] == '+') && - (node.direction[1] == 'x' || node.direction[1] == 'y' || - node.direction[1] == 'z' || node.direction[1] == 'r'), - "Surface current direction is not correctly formatted!"); - mfem::Array attr_marker; mesh::AttrToMarker(h1_fespace.GetParMesh()->bdr_attributes.Max(), node.attributes, attr_marker); @@ -37,7 +30,7 @@ SurfaceCurrentData::SurfaceCurrentData(const config::SurfaceCurrentData &data, else { elems.push_back( - std::make_unique(node.direction, attr_marker, h1_fespace)); + std::make_unique(node.normal, attr_marker, h1_fespace)); } } } diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 2fbf74c6b..8dcc6f06b 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -728,6 +728,108 @@ void ImpedanceBoundaryData::SetUp(json &boundaries) } } +namespace +{ + +// Helper function for extracting a ComponentNode from a json, if the is_port +// value is set to true, will extract the normal vector from either the provided +// keyword argument or from a specified 3 vector. In extracting the normal +// various checks are performed for validity of the input combinations. +auto ParseComponentNode(json &j, bool is_port = true) +{ + ComponentNode node; + node.attributes = j.at("Attributes").get>(); // Required + + try + { + node.direction = j.at("Direction").get(); + } + catch (json::exception) + { + try + { + node.normal = j.at("Direction").get>(); + } + catch (json::exception) + { + MFEM_VERIFY(!is_port, "A port requires the direction be specified"); + return node; + } + } + + for (auto &c : node.direction) + { + c = std::tolower(c); + } + + const auto xpos = node.direction.find("x"); + const auto ypos = node.direction.find("y"); + const auto zpos = node.direction.find("z"); + const auto rpos = node.direction.find("r"); + + const bool xfound = xpos != std::string::npos; + const bool yfound = ypos != std::string::npos; + const bool zfound = zpos != std::string::npos; + const bool rfound = rpos != std::string::npos; + + // Either a keyword direction is specified xor the magnitude is non-zero + MFEM_VERIFY((xfound || yfound || zfound || rfound) ^ (node.NormalMagnitude() > 0.0), + "Keyword and vector specification are mutually exclusive"); + + if (xfound) + { + MFEM_VERIFY(node.direction.length() == 1 || node.direction[xpos - 1] == '-' || + node.direction[xpos - 1] == '+', + "Missing required sign specification on \"X\""); + node.normal[0] = + (node.direction.length() == 1 || node.direction[xpos - 1] == '+') ? 1 : -1; + } + if (yfound) + { + MFEM_VERIFY(node.direction.length() == 1 || node.direction[ypos - 1] == '-' || + node.direction[ypos - 1] == '+', + "Missing sign specification on \"Y\""); + node.normal[1] = + node.direction.length() == 1 || node.direction[ypos - 1] == '+' ? 1 : -1; + } + if (zfound) + { + MFEM_VERIFY(node.direction.length() == 1 || node.direction[zpos - 1] == '-' || + node.direction[zpos - 1] == '+', + "Missing sign specification on \"Z\""); + node.normal[2] = + node.direction.length() == 1 || node.direction[zpos - 1] == '+' ? 1 : -1; + } + if (rfound) + { + if (node.direction.length() == 1) + { + node.direction = "+r"; + } + MFEM_VERIFY(node.direction[rpos - 1] == '-' || node.direction[rpos - 1] == '+', + "Missing sign specification on \"R\""); + MFEM_VERIFY(!xfound && !yfound && !zfound, + "\"R\" cannot be combined with \"X\", \"Y\" or \"Z\""); + } + + double mag = node.NormalMagnitude(); + + double constexpr tol = 1e-6; + if (mag > 0 && std::abs(mag - 1) > tol) + { + std::cout << "Renormalized port direction to: "; + for (auto &x : node.normal) + { + x /= mag; + std::cout << x << ' '; + } + std::cout << '\n'; + } + + return node; +} +} // namespace + void LumpedPortBoundaryData::SetUp(json &boundaries) { auto port = boundaries.find("LumpedPort"); @@ -751,65 +853,51 @@ void LumpedPortBoundaryData::SetUp(json &boundaries) MFEM_VERIFY( port->is_array(), "\"LumpedPort\" and \"Terminal\" should specify an array in the configuration file!"); - for (auto it = port->begin(); it != port->end(); ++it) + for (auto &p : *port) { MFEM_VERIFY( - it->find("Index") != it->end(), + p.find("Index") != p.end(), "Missing \"LumpedPort\" or \"Terminal\" boundary \"Index\" in configuration file!"); - auto ret = mapdata.insert(std::make_pair(it->at("Index"), LumpedPortData())); + auto ret = mapdata.insert(std::make_pair(p.at("Index"), LumpedPortData())); MFEM_VERIFY(ret.second, "Repeated \"Index\" found when processing \"LumpedPort\" or " "\"Terminal\" boundaries in configuration file!"); LumpedPortData &data = ret.first->second; - data.R = it->value("R", data.R); - data.L = it->value("L", data.L); - data.C = it->value("C", data.C); - data.Rs = it->value("Rs", data.Rs); - data.Ls = it->value("Ls", data.Ls); - data.Cs = it->value("Cs", data.Cs); - data.excitation = it->value("Excitation", data.excitation); - if (it->find("Attributes") != it->end()) + data.R = p.value("R", data.R); + data.L = p.value("L", data.L); + data.C = p.value("C", data.C); + data.Rs = p.value("Rs", data.Rs); + data.Ls = p.value("Ls", data.Ls); + data.Cs = p.value("Cs", data.Cs); + data.excitation = p.value("Excitation", data.excitation); + if (p.find("Attributes") != p.end()) { - MFEM_VERIFY(it->find("Elements") == it->end(), + MFEM_VERIFY(p.find("Elements") == p.end(), "Cannot specify both top-level \"Attributes\" list and \"Elements\" for " "\"LumpedPort\" or \"Terminal\" boundary in configuration file!"); - data.nodes.resize(1); - LumpedPortData::Node &node = data.nodes.back(); - node.attributes = it->at("Attributes").get>(); // Required - node.direction = it->value("Direction", node.direction); - if (terminal == boundaries.end()) - { - // Only check direction for true ports, not terminals. - CheckDirection(node.direction, true); - } + + data.nodes.clear(); + data.nodes.emplace_back(ParseComponentNode(p, terminal == boundaries.end())); } else { - auto elements = it->find("Elements"); - MFEM_VERIFY(elements != it->end(), + auto elements = p.find("Elements"); + MFEM_VERIFY(elements != p.end(), "Missing top-level \"Attributes\" list or \"Elements\" for " "\"LumpedPort\" or \"Terminal\" boundary in configuration file!"); - for (auto elem_it = elements->begin(); elem_it != elements->end(); ++elem_it) + for (auto &elem : *elements) { - MFEM_VERIFY(elem_it->find("Attributes") != elem_it->end(), + MFEM_VERIFY(elem.find("Attributes") != elem.end(), "Missing \"Attributes\" list for \"LumpedPort\" or \"Terminal\" " "boundary element in configuration file!"); - data.nodes.emplace_back(); - LumpedPortData::Node &node = data.nodes.back(); - node.attributes = elem_it->at("Attributes").get>(); // Required - node.direction = elem_it->value("Direction", node.direction); - if (terminal == boundaries.end()) - { - // Only check direction for true ports, not terminals. - CheckDirection(node.direction, true); - } + data.nodes.emplace_back(ParseComponentNode(elem, terminal == boundaries.end())); // Cleanup - elem_it->erase("Attributes"); - elem_it->erase("Direction"); - MFEM_VERIFY(elem_it->empty(), + elem.erase("Attributes"); + elem.erase("Direction"); + MFEM_VERIFY(elem.empty(), "Found an unsupported configuration file keyword under \"LumpedPort\" " "or \"Terminal\" boundary element!\n" - << elem_it->dump(2)); + << elem.dump(2)); } } @@ -829,20 +917,20 @@ void LumpedPortBoundaryData::SetUp(json &boundaries) // } // Cleanup - it->erase("Index"); - it->erase("R"); - it->erase("L"); - it->erase("C"); - it->erase("Rs"); - it->erase("Ls"); - it->erase("Cs"); - it->erase("Excitation"); - it->erase("Attributes"); - it->erase("Direction"); - it->erase("Elements"); - MFEM_VERIFY(it->empty(), "Found an unsupported configuration file keyword under " - "\"LumpedPort\" or \"Terminal\"!\n" - << it->dump(2)); + p.erase("Index"); + p.erase("R"); + p.erase("L"); + p.erase("C"); + p.erase("Rs"); + p.erase("Ls"); + p.erase("Cs"); + p.erase("Excitation"); + p.erase("Attributes"); + p.erase("Direction"); + p.erase("Elements"); + MFEM_VERIFY(p.empty(), "Found an unsupported configuration file keyword under " + "\"LumpedPort\" or \"Terminal\"!\n" + << p.dump(2)); } } @@ -901,50 +989,43 @@ void SurfaceCurrentBoundaryData::SetUp(json &boundaries) } MFEM_VERIFY(source->is_array(), "\"SurfaceCurrent\" should specify an array in the configuration file!"); - for (auto it = source->begin(); it != source->end(); ++it) + for (auto &s : *source) { - MFEM_VERIFY(it->find("Index") != it->end(), + MFEM_VERIFY(s.find("Index") != s.end(), "Missing \"SurfaceCurrent\" source \"Index\" in configuration file!"); - auto ret = mapdata.insert(std::make_pair(it->at("Index"), SurfaceCurrentData())); + auto ret = mapdata.insert(std::make_pair(s.at("Index"), SurfaceCurrentData())); MFEM_VERIFY(ret.second, "Repeated \"Index\" found when processing \"SurfaceCurrent\" " "boundaries in configuration file!"); SurfaceCurrentData &data = ret.first->second; - if (it->find("Attributes") != it->end()) + if (s.find("Attributes") != s.end()) { - MFEM_VERIFY(it->find("Elements") == it->end(), + MFEM_VERIFY(s.find("Elements") == s.end(), "Cannot specify both top-level \"Attributes\" list and \"Elements\" for " "\"SurfaceCurrent\" boundary in configuration file!"); - data.nodes.resize(1); - SurfaceCurrentData::Node &node = data.nodes.back(); - node.attributes = it->at("Attributes").get>(); // Required - node.direction = it->value("Direction", node.direction); - CheckDirection(node.direction, true); + data.nodes.clear(); + data.nodes.emplace_back(ParseComponentNode(s)); } else { - auto elements = it->find("Elements"); + auto elements = s.find("Elements"); MFEM_VERIFY( - elements != it->end(), + elements != s.end(), "Missing top-level \"Attributes\" list or \"Elements\" for \"SurfaceCurrent\" " "boundary in configuration file!"); - for (auto elem_it = elements->begin(); elem_it != elements->end(); ++elem_it) + for (auto elem : *elements) { MFEM_VERIFY( - elem_it->find("Attributes") != elem_it->end(), + elem.find("Attributes") != elem.end(), "Missing \"Attributes\" list for \"SurfaceCurrent\" boundary element in " "configuration file!"); - data.nodes.emplace_back(); - SurfaceCurrentData::Node &node = data.nodes.back(); - node.attributes = it->at("Attributes").get>(); // Required - node.direction = it->value("Direction", node.direction); - CheckDirection(node.direction, true); + data.nodes.emplace_back(ParseComponentNode(s)); // Cleanup - elem_it->erase("Attributes"); - elem_it->erase("Direction"); - MFEM_VERIFY(elem_it->empty(), "Found an unsupported configuration file keyword " - "under \"SurfaceCurrent\" boundary element!\n" - << elem_it->dump(2)); + elem.erase("Attributes"); + elem.erase("Direction"); + MFEM_VERIFY(elem.empty(), "Found an unsupported configuration file keyword " + "under \"SurfaceCurrent\" boundary element!\n" + << elem.dump(2)); } } @@ -957,14 +1038,14 @@ void SurfaceCurrentBoundaryData::SetUp(json &boundaries) // } // Cleanup - it->erase("Index"); - it->erase("Attributes"); - it->erase("Direction"); - it->erase("Elements"); + s.erase("Index"); + s.erase("Attributes"); + s.erase("Direction"); + s.erase("Elements"); MFEM_VERIFY( - it->empty(), + s.empty(), "Found an unsupported configuration file keyword under \"SurfaceCurrent\"!\n" - << it->dump(2)); + << s.dump(2)); } } diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index b3e785565..166f571d4 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -4,6 +4,7 @@ #ifndef PALACE_UTILS_CONFIG_FILE_HPP #define PALACE_UTILS_CONFIG_FILE_HPP +#include #include #include #include @@ -345,6 +346,40 @@ struct ImpedanceBoundaryData : public internal::DataVector void SetUp(json &boundaries); }; +// A component Node consists of a list of attributes making up a single +// element of a potentially multielement node, and a direction and/or a normal +// defining the incident field. These are used for Lumped Ports, Terminals, +// and Surface Currents. +// +// Exactly one of either the direction or the normal direction must be defined. +struct ComponentNode +{ + // Optional string defining source excitation field direction. Options are "X", "Y", "Z", + // or "R", preceeded with a "+" or "-" for the direction. The first three options are for + // uniform lumped ports and the last is for a coaxial lumped port (radial excitation). + std::string direction = ""; + + // Vector defining the normal direction for this port. Used in Cartesian + // coordinate system ports, with "X", "Y", and "Z" mapping to (1,0,0), + // (0,1,0), and (0,0,1) respectively. Any user specified value is normalized + // to a tolerance of 1e-6, and if normalization is needed, printed to the terminal. + std::array normal{{0.0, 0.0, 0.0}}; + + // List of boundary attributes for this lumped port element. + std::vector attributes = {}; + + // Convenience function for computing the normal magnitude. + inline double NormalMagnitude() const + { + double mag = 0.0; + for (auto x : normal) + { + mag += x * x; + } + return std::sqrt(mag); + } +}; + struct LumpedPortData { public: @@ -361,19 +396,7 @@ struct LumpedPortData // Flag for source term in driven and transient simulations. bool excitation = false; - // For each lumped port index, each Node contains a list of attributes making up a single - // element of a potentially multielement port. - struct Node - { - // String defining source excitation field direction. Options are "X", "Y", "Z", or "R", - // preceeded with a "+" or "-" for the direction. The first three options are for - // uniform lumped ports and the last is for a coaxial lumped port (radial excitation). - std::string direction = ""; - - // List of boundary attributes for this lumped port element. - std::vector attributes = {}; - }; - std::vector nodes = {}; + std::vector nodes = {}; }; struct LumpedPortBoundaryData : public internal::DataMap @@ -409,17 +432,7 @@ struct SurfaceCurrentData public: // For each surface current source index, each Node contains a list of attributes making // up a single element of a potentially multielement current source. - struct Node - { - // String defining surface current source excitation direction. Options are "X", "Y", - // "Z", or "R", preceeded with a "+" or "-" for the direction. The first three options - // are for uniform sources and the last is for a coaxial source (radial excitation). - std::string direction = ""; - - // List of boundary attributes for this surface current element. - std::vector attributes = {}; - }; - std::vector nodes = {}; + std::vector nodes = {}; }; struct SurfaceCurrentBoundaryData : public internal::DataMap diff --git a/scripts/schema/config/boundaries.json b/scripts/schema/config/boundaries.json index 442edf303..2ada4e0a5 100644 --- a/scripts/schema/config/boundaries.json +++ b/scripts/schema/config/boundaries.json @@ -117,7 +117,7 @@ "additionalItems": false, "items": { "type": "integer"} }, - "Direction": { "type": "string" }, + "Direction": { "$ref": "#/$defs/Direction" }, "R": { "type": "number" }, "L": { "type": "number" }, "C": { "type": "number" }, @@ -142,7 +142,7 @@ "additionalItems": false, "items": { "type": "integer"} }, - "Direction": { "type": "string" } + "Direction": { "$ref": "#/$defs/Direction" } } } } @@ -206,7 +206,7 @@ "additionalItems": false, "items": { "type": "integer"} }, - "Direction": { "type": "string" }, + "Direction": { "$ref": "#/$defs/Direction" }, "Elements": { "type": "array", @@ -224,7 +224,7 @@ "additionalItems": false, "items": { "type": "integer"} }, - "Direction": { "type": "string" } + "Direction": { "$ref": "#/$defs/Direction" } } } } @@ -328,7 +328,7 @@ "additionalItems": false, "items": { "type": "integer"} }, - "Direction": { "type": "string", "minLength": 1, "maxLength": 2 } + "Direction": { "$ref": "#/$defs/Direction" } } } }, @@ -383,5 +383,25 @@ } } } + }, + "$defs": + { + "Direction": + { + "anyOf": + [ + { + "type": "string", + "enum": ["R", "X", "Y", "Z", "+R", "+X", "+Y", "+Z", "-R", "-X", "-Y", "-Z", + "r", "x", "y", "z", "+r", "+x", "+y", "+z", "-r", "-x", "-y", "-z"] + }, + { + "type": "array", + "items": { "type": "number" }, + "minItems": 3, "maxItems": 3 + } + ] + } } } +