Skip to content

Commit 46e4f85

Browse files
laylagilaylagi
committed
Initial commit for adding Electrical Current Dipole
1 parent 8a5367f commit 46e4f85

File tree

11 files changed

+452
-17
lines changed

11 files changed

+452
-17
lines changed

‎examples/cylinder/driven_dipole.json‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
{
2+
"Problem":
3+
{
4+
"Type": "Driven",
5+
"Verbose": 2,
6+
"Output": "postpro/driven_dipole"
7+
},
8+
"Model":
9+
{
10+
"Mesh": "mesh/cylinder_hex.msh",
11+
"L0": 1.0e-2
12+
},
13+
"Domains":
14+
{
15+
"Materials":
16+
[
17+
{
18+
"Attributes": [1],
19+
"Permeability": 1.0,
20+
"Permittivity": 2.08,
21+
"LossTan": 0.0004
22+
}
23+
],
24+
"CurrentDipole": [
25+
{
26+
"Index": 1,
27+
"Moment": [2.0, 0.0, 0.0],
28+
"Center": [0.0, 0.0, 2.0]
29+
},
30+
{
31+
"Index": 2,
32+
"Moment": [0.0, -2.0, 0.0],
33+
"Center": [0.0, 0.0, 3.48]
34+
}
35+
],
36+
"Postprocessing":
37+
{
38+
"Energy":
39+
[
40+
{
41+
"Index": 1,
42+
"Attributes": [1]
43+
}
44+
]
45+
}
46+
},
47+
"Boundaries":
48+
{
49+
"PEC":
50+
{
51+
"Attributes": [2, 3, 4]
52+
}
53+
},
54+
"Solver":
55+
{
56+
"Order": 4,
57+
"Device": "CPU",
58+
"Driven":
59+
{
60+
"MinFreq": 1.0,
61+
"MaxFreq": 4.0,
62+
"FreqStep": 0.25,
63+
"SaveStep": 2
64+
},
65+
"Linear":
66+
{
67+
"Type": "Default",
68+
"KSPType": "GMRES",
69+
"Tol": 1.0e-9,
70+
"MaxIts": 100
71+
}
72+
}
73+
}

‎palace/models/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
target_sources(${LIB_TARGET_NAME}
99
PRIVATE
1010
${CMAKE_CURRENT_SOURCE_DIR}/curlcurloperator.cpp
11+
${CMAKE_CURRENT_SOURCE_DIR}/currentdipoleoperator.cpp
1112
${CMAKE_CURRENT_SOURCE_DIR}/domainpostoperator.cpp
1213
${CMAKE_CURRENT_SOURCE_DIR}/farfieldboundaryoperator.cpp
1314
${CMAKE_CURRENT_SOURCE_DIR}/laplaceoperator.cpp

‎palace/models/currentdipoleoperator.cpp‎

Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
// Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2+
// SPDX-License-Identifier: Apache-2.0
3+
4+
#include "currentdipoleoperator.hpp"
5+
6+
#include "fem/coefficient.hpp"
7+
#include "utils/communication.hpp"
8+
#include "utils/geodata.hpp"
9+
#include "utils/iodata.hpp"
10+
11+
namespace palace
12+
{
13+
14+
CurrentDipoleData::CurrentDipoleData(const config::CurrentDipoleData &data,
15+
const mfem::ParMesh &mesh)
16+
{
17+
// Set up dipole moment vector
18+
moment.SetSize(mesh.SpaceDimension());
19+
MFEM_VERIFY(data.moment.size() == static_cast<std::size_t>(mesh.SpaceDimension()),
20+
"Current dipole moment dimension must match mesh space dimension!");
21+
for (int d = 0; d < mesh.SpaceDimension(); d++)
22+
{
23+
moment[d] = data.moment[d];
24+
}
25+
26+
// Set up dipole center position
27+
center.SetSize(mesh.SpaceDimension());
28+
MFEM_VERIFY(data.center.size() == static_cast<std::size_t>(mesh.SpaceDimension()),
29+
"Current dipole center dimension must match mesh space dimension!");
30+
for (int d = 0; d < mesh.SpaceDimension(); d++)
31+
{
32+
center[d] = data.center[d];
33+
}
34+
35+
// Create the VectorDeltaCoefficient
36+
dipole_coeff = std::make_unique<mfem::VectorDeltaCoefficient>(moment);
37+
dipole_coeff->SetDeltaCenter(center);
38+
}
39+
40+
mfem::VectorDeltaCoefficient *CurrentDipoleData::GetDeltaCoefficient(double scale) const
41+
{
42+
// Scale the dipole coefficient by the given factor
43+
dipole_coeff->SetScale(scale);
44+
return dipole_coeff.get();
45+
}
46+
47+
CurrentDipoleOperator::CurrentDipoleOperator(const IoData &iodata,
48+
const mfem::ParMesh &mesh)
49+
{
50+
// Set up current dipole source properties.
51+
SetUpDipoleProperties(iodata, mesh);
52+
PrintDipoleInfo(iodata, mesh);
53+
}
54+
55+
void CurrentDipoleOperator::SetUpDipoleProperties(const IoData &iodata,
56+
const mfem::ParMesh &mesh)
57+
{
58+
// Set up current dipole data structures.
59+
for (const auto &[idx, data] : iodata.domains.current_dipole)
60+
{
61+
dipoles.emplace(idx, CurrentDipoleData(data, mesh));
62+
}
63+
}
64+
65+
void CurrentDipoleOperator::PrintDipoleInfo(const IoData &iodata, const mfem::ParMesh &mesh)
66+
{
67+
if (dipoles.empty())
68+
{
69+
return;
70+
}
71+
72+
Mpi::Print("\nConfiguring electrical current dipole source terms:\n");
73+
for (const auto &[idx, data] : dipoles)
74+
{
75+
// Convert center coordinates back to physical units for display
76+
mfem::Vector physical_center = data.center;
77+
iodata.units.DimensionalizeInPlace<Units::ValueType::LENGTH>(physical_center);
78+
Mpi::Print(" Dipole {:d}: \n"
79+
" \tCenter = ({:.3e}, {:.3e}, {:.3e})\n"
80+
" \tMoment = ({:.3e}, {:.3e}, {:.3e})\n",
81+
idx, physical_center[0], physical_center[1],
82+
(mesh.SpaceDimension() > 2) ? physical_center[2] : 0.0, data.moment[0],
83+
data.moment[1], (mesh.SpaceDimension() > 2) ? data.moment[2] : 0.0);
84+
}
85+
}
86+
87+
const CurrentDipoleData &CurrentDipoleOperator::GetDipole(int idx) const
88+
{
89+
return dipoles.at(idx);
90+
}
91+
92+
std::vector<mfem::Vector> CurrentDipoleOperator::GetMoments() const
93+
{
94+
std::vector<mfem::Vector> moments;
95+
moments.reserve(dipoles.size());
96+
for (const auto &[idx, dipole] : dipoles)
97+
{
98+
moments.push_back(dipole.moment);
99+
}
100+
return moments;
101+
}
102+
103+
std::vector<mfem::Vector> CurrentDipoleOperator::GetCenters() const
104+
{
105+
std::vector<mfem::Vector> centers;
106+
centers.reserve(dipoles.size());
107+
for (const auto &[idx, dipole] : dipoles)
108+
{
109+
centers.push_back(dipole.center);
110+
}
111+
return centers;
112+
}
113+
114+
void CurrentDipoleOperator::AddExcitationDomainIntegrators(int idx, mfem::LinearForm &rhs)
115+
{
116+
// Add only the specific dipole with the given index
117+
auto it = dipoles.find(idx);
118+
if (it != dipoles.end())
119+
{
120+
const auto &dipole = it->second;
121+
122+
// Create a VectorDeltaCoefficient for this specific dipole (RHS = -J_source)
123+
auto dipole_coeff = std::make_unique<mfem::VectorDeltaCoefficient>(dipole.moment);
124+
dipole_coeff->SetDeltaCenter(dipole.center);
125+
dipole_coeff->SetScale(-1.0);
126+
127+
// Add as domain integrator
128+
rhs.AddDomainIntegrator(new mfem::VectorFEDomainLFIntegrator(*dipole_coeff));
129+
130+
// Store the coefficient to prevent deletion.
131+
dipole_integrator_coeffs.push_back(std::move(dipole_coeff));
132+
}
133+
}
134+
135+
void CurrentDipoleOperator::AddExcitationDomainIntegrators(mfem::LinearForm &rhs)
136+
{
137+
// Add each dipole as a separate integrator
138+
for (const auto &[idx, dipole] : dipoles)
139+
{
140+
AddExcitationDomainIntegrators(idx, rhs);
141+
}
142+
}
143+
144+
} // namespace palace

‎palace/models/currentdipoleoperator.hpp‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
// Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2+
// SPDX-License-Identifier: Apache-2.0
3+
4+
#ifndef PALACE_MODELS_CURRENT_DIPOLE_OPERATOR_HPP
5+
#define PALACE_MODELS_CURRENT_DIPOLE_OPERATOR_HPP
6+
7+
#include <map>
8+
#include <memory>
9+
#include <vector>
10+
#include <mfem.hpp>
11+
12+
namespace palace
13+
{
14+
15+
class IoData;
16+
class SumVectorCoefficient;
17+
18+
namespace config
19+
{
20+
21+
struct CurrentDipoleData;
22+
23+
} // namespace config
24+
25+
// Helper class for current dipole sources in a model.
26+
class CurrentDipoleData
27+
{
28+
public:
29+
// Current dipole properties.
30+
mfem::Vector moment; // Dipole moment vector
31+
mfem::Vector center; // Dipole center position
32+
33+
// Internal MFEM coefficient for the dipole
34+
std::unique_ptr<mfem::VectorDeltaCoefficient> dipole_coeff;
35+
36+
public:
37+
CurrentDipoleData(const config::CurrentDipoleData &data, const mfem::ParMesh &mesh);
38+
39+
// Get the VectorDeltaCoefficient for this dipole
40+
// TODO: scale defaults to 1. Do we need to expose this?
41+
mfem::VectorDeltaCoefficient *GetDeltaCoefficient(double scale = 1.0) const;
42+
};
43+
44+
//
45+
// A class handling current dipole sources and their excitation.
46+
//
47+
class CurrentDipoleOperator
48+
{
49+
private:
50+
// Mapping from dipole index to data structure containing dipole information.
51+
std::map<int, CurrentDipoleData> dipoles;
52+
53+
// Storage for integrator coefficients to manage their lifetime
54+
mutable std::vector<std::unique_ptr<mfem::VectorDeltaCoefficient>>
55+
dipole_integrator_coeffs;
56+
57+
void SetUpDipoleProperties(const IoData &iodata, const mfem::ParMesh &mesh);
58+
void PrintDipoleInfo(const IoData &iodata, const mfem::ParMesh &mesh);
59+
60+
public:
61+
CurrentDipoleOperator(const IoData &iodata, const mfem::ParMesh &mesh);
62+
63+
// Access data structures for current dipoles.
64+
const CurrentDipoleData &GetDipole(int idx) const;
65+
auto begin() const { return dipoles.begin(); }
66+
auto end() const { return dipoles.end(); }
67+
auto rbegin() const { return dipoles.rbegin(); }
68+
auto rend() const { return dipoles.rend(); }
69+
auto Size() const { return dipoles.size(); }
70+
bool Empty() const { return dipoles.empty(); }
71+
72+
// Get moment and center arrays separately
73+
std::vector<mfem::Vector> GetMoments() const;
74+
std::vector<mfem::Vector> GetCenters() const;
75+
76+
// Add integrators directly to a LinearForm
77+
void AddExcitationDomainIntegrators(mfem::LinearForm &rhs);
78+
void AddExcitationDomainIntegrators(int idx, mfem::LinearForm &rhs);
79+
};
80+
81+
} // namespace palace
82+
83+
#endif // PALACE_MODELS_CURRENT_DIPOLE_OPERATOR_HPP

‎palace/models/portexcitations.cpp‎

Lines changed: 29 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33

44
#include "portexcitations.hpp"
55

6+
#include "currentdipoleoperator.hpp"
67
#include "lumpedportoperator.hpp"
78
#include "surfacecurrentoperator.hpp"
89
#include "waveportoperator.hpp"
@@ -40,6 +41,11 @@ namespace palace
4041
to(" Surface current port{} {:2d}\n", (ex.current_port.size() > 1) ? "s" : "",
4142
fmt::join(ex.current_port, " "));
4243
}
44+
if (!ex.current_dipole.empty())
45+
{
46+
to(" Current dipole{} {:2d}\n", (ex.current_dipole.size() > 1) ? "s" : "",
47+
fmt::join(ex.current_dipole, " "));
48+
}
4349
i++;
4450
}
4551
return fmt::to_string(buf);
@@ -49,14 +55,16 @@ void to_json(nlohmann::json &j, const PortExcitations::SingleExcitationSpec &p)
4955
{
5056
j = nlohmann::json{{"LumpedPort", p.lumped_port},
5157
{"WavePort", p.wave_port},
52-
{"SurfaceCurrent", p.current_port}};
58+
{"SurfaceCurrent", p.current_port},
59+
{"CurrentDipole", p.current_dipole}};
5360
}
5461

5562
void from_json(const nlohmann::json &j, PortExcitations::SingleExcitationSpec &p)
5663
{
5764
j.at("LumpedPort").get_to(p.lumped_port);
5865
j.at("WavePort").get_to(p.wave_port);
5966
j.at("SurfaceCurrent").get_to(p.current_port);
67+
j.at("CurrentDipole").get_to(p.current_dipole);
6068
}
6169

6270
void to_json(nlohmann::json &j, const PortExcitations &p)
@@ -71,7 +79,8 @@ void from_json(const nlohmann::json &j, PortExcitations &p)
7179

7280
PortExcitations::PortExcitations(const LumpedPortOperator &lumped_port_op,
7381
const WavePortOperator &wave_port_op,
74-
const SurfaceCurrentOperator &surf_j_op)
82+
const SurfaceCurrentOperator &surf_j_op,
83+
const CurrentDipoleOperator &dipole_op)
7584
{
7685
for (const auto &[idx, port] : lumped_port_op)
7786
{
@@ -106,6 +115,24 @@ PortExcitations::PortExcitations(const LumpedPortOperator &lumped_port_op,
106115
ex_spec.current_port = current_port_idx;
107116
}
108117
}
118+
119+
// Current dipoles are always excited. Add them to all existing excitations.
120+
std::vector<int> current_dipole_idx;
121+
for (const auto &[idx, dipole] : dipole_op)
122+
{
123+
current_dipole_idx.push_back(idx);
124+
}
125+
if (!current_dipole_idx.empty())
126+
{
127+
if (excitations.empty())
128+
{
129+
excitations.try_emplace(1, SingleExcitationSpec{});
130+
}
131+
for (auto &[ex_idx, ex_spec] : excitations)
132+
{
133+
ex_spec.current_dipole = current_dipole_idx;
134+
}
135+
}
109136
};
110137

111138
} // namespace palace

0 commit comments

Comments
 (0)