Mesh rescaling finished, line integration part 1.

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

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;
 };
-

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);
+}

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,25 +194,18 @@ int main(int argc, const char *argv[])
     gmsh::open( mpl.sim_gmshmodel() );
 
     model mod( mpl.sim_geomorder(), mpl.sim_approxorder() );
-    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 };
+        mod.build(sf);
+    else
+        mod.build();
 
-        gmsh::model::mesh::affineTransform(tr);
+    if ( not mpl.validate_model_params(mod) )
+    {
+        std::cout << "Configuration problem, exiting" << std::endl;
+        return 1;
     }
-    */
 
     if ( mpl.sim_usegpu() )
     {

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();

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;

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;

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;

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;

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;

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();