diff --git a/doc/lua_api.md b/doc/lua_api.md
index d74839f0d9b10cdf2905ff889d9b0c4ee0cb0a05..67c03cca57ed28f74775bc29e649c3cf5fbcab86 100644
--- a/doc/lua_api.md
+++ b/doc/lua_api.md
@@ -27,6 +27,9 @@ This document describes the API available on the version v0.1 of the solver.
- `postpro.silo_output_rate` (integer):
- `postpro.cycle_print_rate` (integer):
+### Mesh variables
+- `mesh.scalefactor` (real): apply a scaling factor to the mesh.
+
### Callable functions
- `enable_boundary_sources()`: takes a boolean parameter that specifies if boundary sources should be enabled or not
- `enable_interface_sources()`: takes a boolean parameter that specifies if the sources applied on internal interfaces should be enabled or not
diff --git a/include/gmsh_io.h b/include/gmsh_io.h
index 7585f304322d8c578dc530c39c87aa35b9b8d060..0824294ad5f41df66b74ac9dd2fd91aacb61348b 100644
--- a/include/gmsh_io.h
+++ b/include/gmsh_io.h
@@ -66,10 +66,15 @@ class model
std::vector<boundary_descriptor> bnd_descriptors;
element_connectivity<element_key> conn;
+ using entofs_pair = std::pair<size_t, size_t>;
+ std::map<size_t, entofs_pair> etag_to_entity_offset;
+
void map_boundaries(void);
void import_gmsh_entities(void);
void update_connectivity(const entity&, size_t);
+ void build_priv(void);
+
public:
model();
model(int, int);
@@ -77,8 +82,8 @@ public:
model(const char *, int, int);
~model();
- void remesh(void);
- void remesh(int);
+ void build(void);
+ void build(double);
const element_connectivity<element_key>& connectivity() const
{
@@ -109,5 +114,5 @@ public:
entity& entity_at(size_t);
const entity& entity_at(size_t) const;
+ entofs_pair lookup_tag(size_t tag) const;
};
-
diff --git a/src/gmsh_io.cpp b/src/gmsh_io.cpp
index 5e2e23cb5ef29783f0204c7ffdbdeb78026d3170..bdddb6e911c5067c8cb38389327dacfd209962a6 100644
--- a/src/gmsh_io.cpp
+++ b/src/gmsh_io.cpp
@@ -35,24 +35,19 @@ basis_grad_name(int order)
model::model()
: geometric_order(1), approximation_order(1)
-{
- remesh(geometric_order);
-}
+{}
model::model(int pgo, int pao)
: geometric_order(pgo), approximation_order(pao)
{
assert(geometric_order > 0);
assert(approximation_order > 0);
-
- remesh(geometric_order);
}
model::model(const char *filename)
: geometric_order(1), approximation_order(1)
{
gmsh::open( std::string(filename) );
- remesh(geometric_order);
}
model::model(const char *filename, int pgo, int pao)
@@ -62,7 +57,6 @@ model::model(const char *filename, int pgo, int pao)
assert(approximation_order > 0);
gmsh::open( std::string(filename) );
- remesh(geometric_order);
}
model::~model()
@@ -71,9 +65,8 @@ model::~model()
}
void
-model::remesh()
+model::build_priv()
{
- gmm::generate( DIMENSION(3) );
gmm::setOrder( geometric_order );
gvp_t gmsh_entities;
@@ -103,10 +96,22 @@ model::remesh()
}
void
-model::remesh(int pgo)
+model::build()
{
- geometric_order = pgo;
- remesh();
+ gmm::generate( DIMENSION(3) );
+ build_priv();
+}
+
+void
+model::build(double sf)
+{
+ gmm::generate( DIMENSION(3) );
+
+ std::vector<double> tr = {sf, 0, 0, 0,
+ 0, sf, 0, 0,
+ 0, 0, sf, 0 };
+ gmm::affineTransform(tr);
+ build_priv();
}
std::vector<entity>::const_iterator
@@ -231,6 +236,17 @@ model::import_gmsh_entities(void)
.gorder = geometric_order, .aorder = approximation_order});
e.sort_by_orientation();
+
+ for (size_t i = 0; i < e.num_cells(); i++)
+ {
+ const auto& pe = e.cell(i);
+ auto eitor = etag_to_entity_offset.find( pe.element_tag() );
+ if ( eitor != etag_to_entity_offset.end() )
+ throw std::logic_error("Duplicate tag");
+
+ etag_to_entity_offset[ pe.element_tag() ] = std::make_pair(entity_index, i);
+ }
+
e.base(dof_base, flux_base, index_base);
e.number(entity_index);
dof_base += e.num_dofs();
@@ -313,3 +329,9 @@ model::map_boundaries(void)
std::cout << normal << " " << interface << " " << boundary << std::endl;
*/
}
+
+model::entofs_pair
+model::lookup_tag(size_t tag) const
+{
+ return etag_to_entity_offset.at(tag);
+}
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index c40793b34c72d20372709c9c8cc2e1699bb8e135..aee2197fbfc8f0e985ac98199b8c36f5fb58ae57 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -66,6 +66,54 @@ void initialize_solver(const model& mod, State& state, const maxwell::parameter_
init_matparams(mod, state, mpl);
}
+class line_integrator
+{
+ std::vector<size_t> element_tags;
+ std::vector<point_3d> sampling_points_phys;
+ std::vector<point_3d> sampling_points_ref;
+ point_3d increment;
+
+public:
+ line_integrator(const point_3d& start,
+ const point_3d& end, size_t samples)
+ {
+ increment = (end - start)/samples;
+ for (size_t i = 0; i < samples; i++)
+ {
+ point_3d sp = start + increment*(i+0.5);
+ sampling_points_phys.push_back(sp);
+
+ size_t etag;
+ int etype;
+ std::vector<size_t> ntags;
+ double u, v, w;
+ gmsh::model::mesh::getElementByCoordinates(sp.x(), sp.y(), sp.z(),
+ etag, etype, ntags, u, v, w, 3, true);
+
+ element_tags.push_back(etag);
+ sampling_points_ref.push_back( point_3d(u,v,w) );
+ }
+ }
+
+ double integrate(const model& mod, const vecxd& dofs)
+ {
+ double integral_val = 0.0;
+ for (size_t i = 0; i < sampling_points_phys.size(); i++)
+ {
+ auto [entnum, ofs] = mod.lookup_tag(element_tags.at(i));
+ auto& e = mod.entity_at(entnum);
+ auto& pe = e.cell(ofs);
+ auto& re = e.cell_refelem(pe);
+ auto dof_ofs = e.cell_global_dof_offset(ofs);
+ auto dof_num = re.num_basis_functions();
+ vecxd phi = re.basis_functions( sampling_points_ref[i] );
+ vecxd sol = dofs.segment(dof_ofs, dof_num);
+ integral_val += sol.dot(phi) * increment.x();
+ }
+
+ return integral_val;
+ }
+};
template<typename State>
void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
@@ -87,6 +135,10 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
std::cout << "s each." << std::endl;
std::cout << "Time integrator: " << mpl.sim_timeIntegratorName() << std::endl;
+ point_3d start(-46e-9, 0.0, 38e-9);
+ point_3d end(46e-9, 0.0, 38e-9);
+ line_integrator lint(start, end, 50);
+
timecounter tc;
tc.tic();
prepare_sources(mod, state, mpl);
@@ -108,6 +160,7 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
std::cout << "Cycle " << i << ": t = " << state.curr_time << " s";
std::cout << ", DOFs/s: " << (mod.num_dofs()*6*cycle_print_rate)/(time) << std::endl;
tc.tic();
+ std::cout << "Voltage = " << lint.integrate(mod, state.emf_curr.Ex) << std::endl;
}
}
}
@@ -141,26 +194,19 @@ int main(int argc, const char *argv[])
gmsh::open( mpl.sim_gmshmodel() );
model mod( mpl.sim_geomorder(), mpl.sim_approxorder() );
+
+ auto [scaled, sf] = mpl.mesh_scalefactor();
+ if (scaled)
+ mod.build(sf);
+ else
+ mod.build();
+
if ( not mpl.validate_model_params(mod) )
{
std::cout << "Configuration problem, exiting" << std::endl;
return 1;
}
- /*
- //must be after gmm::generate() and before gmm::setOrder() :(
-
- auto [scaled, sf] = mpl.mesh_scalefactor();
- if (scaled)
- {
- std::vector<double> tr = {sf, 0, 0, 0,
- 0, sf, 0, 0,
- 0, 0, sf, 0 };
-
- gmsh::model::mesh::affineTransform(tr);
- }
- */
-
if ( mpl.sim_usegpu() )
{
#ifdef ENABLE_GPU_SOLVER
diff --git a/tests/profile_differentiation_gpu.cpp b/tests/profile_differentiation_gpu.cpp
index 132650196d09965e0bb3d47a7f59a8528a1122a0..501a8af1894c3afa284f345c2fedf0ca6b453460 100644
--- a/tests/profile_differentiation_gpu.cpp
+++ b/tests/profile_differentiation_gpu.cpp
@@ -48,6 +48,7 @@ int profile_differentiation(int geometric_order, int approximation_order)
double h = 0.04;
make_geometry(0,h);
model mod(geometric_order, approximation_order);
+ mod.build();
auto& e0 = mod.entity_at(0);
e0.sort_by_orientation();
diff --git a/tests/test_basics.cpp b/tests/test_basics.cpp
index 5afd04561f9e86f8c0c60c1f4c0a25b151e5c156..3b3a0e9b58d41585abb600ccf88c01870a44f584 100644
--- a/tests/test_basics.cpp
+++ b/tests/test_basics.cpp
@@ -41,6 +41,7 @@ test_basics(int geometric_order, int approximation_order)
make_geometry(0, 0.1);
model mod(geometric_order, approximation_order);
+ mod.build();
auto idx = geometric_order - 1;
@@ -163,6 +164,7 @@ test_normals(int geometric_order, int approximation_order)
{
make_geometry(0, 0.1);
model mod(geometric_order, approximation_order);
+ mod.build();
gmsh::vectorpair pgs;
gm::getPhysicalGroups(pgs);
@@ -246,6 +248,7 @@ test_normals_2()
{
make_geometry(0, 0.1);
model mod(1, 1);
+ mod.build();
vec3d totflux = vec3d::Zero();
double tf = 0.0;
diff --git a/tests/test_differentiation.cpp b/tests/test_differentiation.cpp
index 4bb0b36f78252f42c289a7db728e68704c621434..884d8f56f95c364f34e08b780a0abf3302a85c86 100644
--- a/tests/test_differentiation.cpp
+++ b/tests/test_differentiation.cpp
@@ -70,6 +70,7 @@ int test_differentiation_convergence(int geometric_order, int approximation_orde
double h = sizes[i];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
+ mod.build();
std::stringstream ss;
ss << "diff_go_" << geometric_order << "_ao_" << approximation_order;
diff --git a/tests/test_differentiation_gpu.cpp b/tests/test_differentiation_gpu.cpp
index 3f00c3d802f2cd5f64ab2d0f1f755abb5f83aa1f..c9afbd298d2a3cf870712fb8bfb71544efaf74a9 100644
--- a/tests/test_differentiation_gpu.cpp
+++ b/tests/test_differentiation_gpu.cpp
@@ -72,6 +72,7 @@ int test_differentiation_convergence(int geometric_order, int approximation_orde
double h = sizes[i];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
+ mod.build();
#ifdef WRITE_TEST_OUTPUTS
std::stringstream ss;
diff --git a/tests/test_lifting.cpp b/tests/test_lifting.cpp
index 3cb36f25d58a44b14146f207e32e1171d56de8b9..a543f9080da4d427ec9c61a3f47678f526bb281a 100644
--- a/tests/test_lifting.cpp
+++ b/tests/test_lifting.cpp
@@ -82,6 +82,7 @@ int test_lifting(int geometric_order, int approximation_order)
double h = sizes[sz];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
+ mod.build();
#ifdef WRITE_TEST_OUTPUTS
std::stringstream ss;
diff --git a/tests/test_lifting_gpu.cpp b/tests/test_lifting_gpu.cpp
index cf8fe0ba9559612ae996fdec132813a8e67ac3ff..cd289c8390a77cf7cc933f8672bde94d7861e48e 100644
--- a/tests/test_lifting_gpu.cpp
+++ b/tests/test_lifting_gpu.cpp
@@ -83,6 +83,7 @@ int test_lifting(int geometric_order, int approximation_order)
double h = sizes[sz];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
+ mod.build();
#ifdef WRITE_TEST_OUTPUTS
std::stringstream ss;
diff --git a/tests/test_mass.cpp b/tests/test_mass.cpp
index fcd49b79ec32dc24909e0f062d9ad1dc71721104..9e93b8f55ba5057e1bb5458d990d2ccaeadf7f54 100644
--- a/tests/test_mass.cpp
+++ b/tests/test_mass.cpp
@@ -40,6 +40,7 @@ int test_mass_convergence(int geometric_order, int approximation_order)
double h = sizes[i];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
+ mod.build();
auto& e0 = mod.entity_at(0);
e0.sort_by_orientation();