Skip to content
Snippets Groups Projects
Commit cdcc53a1 authored by Matteo Cicuttin's avatar Matteo Cicuttin
Browse files

Got rid of dof_range, kept track of model/world offsets directly in entity.

parent 6359a404
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 207 additions and 164 deletions
include/entity.h
+ 21
− 8
View file @ cdcc53a1
......@@ -27,10 +27,18 @@ class entity
int parent_dim;
int parent_tag;
size_t m_dof_base;
size_t m_flux_base;
size_t m_index_base;
size_t m_entity_number;
/* Each entity has model offsets, which are offsets relative to the
* computational unit, and world offsets, which are offsets needed
* to recover the whole-world vectors of DOFs. */
size_t m_dof_base_model;
size_t m_flux_base_model;
size_t m_index_base_model;
size_t m_entity_number_model;
size_t m_dof_base_world;
size_t m_flux_base_world;
size_t m_index_base_world;
size_t m_entity_number_world;
entity_ordering cur_elem_ordering;
......@@ -82,11 +90,14 @@ public:
const physical_element& cell(size_t) const;
const reference_element& cell_refelem(size_t) const;
const reference_element& cell_refelem(const physical_element&) const;
size_t cell_global_index(size_t) const;
size_t cell_model_index(size_t) const;
size_t cell_world_index(size_t) const;
size_t cell_local_index_by_gmsh(size_t) const;
size_t cell_global_index_by_gmsh(size_t) const;
size_t cell_model_index_by_gmsh(size_t) const;
size_t cell_world_index_by_gmsh(size_t) const;
size_t cell_local_dof_offset(size_t) const;
size_t cell_global_dof_offset(size_t) const;
size_t cell_model_dof_offset(size_t) const;
size_t cell_world_dof_offset(size_t) const;
size_t num_faces(void) const;
std::vector<size_t> face_tags() const;
......@@ -96,8 +107,10 @@ public:
size_t num_dofs(void) const;
size_t num_fluxes(void) const;
void base(size_t, size_t, size_t);
void base_model(size_t, size_t, size_t);
void base_world(size_t, size_t, size_t);
size_t dof_base(void) const;
size_t dof_base_world(void) const;
size_t flux_base(void) const;
size_t index_base(void) const;
size_t number(void) const;
......
include/gmsh_io.h
+ 30
− 44
View file @ cdcc53a1
......@@ -82,42 +82,23 @@ struct boundary_descriptor
#define IMPLIES(a, b) ((not(a)) or (b))
#define IFF(a,b) (IMPLIES(a,b) and IMPLIES(b,a))
#ifdef USE_MPI
struct dof_range
{
int tag;
int dim;
size_t base;
std::vector<size_t> dof_offsets;
//dof_range(dof_range&&) = default;
dof_range(const entity& e, size_t b)
{
tag = e.gmsh_tag();
dim = e.gmsh_dim();
base = b;
dof_offsets.resize( e.num_dofs() );
auto dofs_per_cell = e.num_dofs()/e.num_cells();
assert( dofs_per_cell*e.num_cells() == e.num_dofs() );
for (size_t iT = 0; iT < e.num_cells(); iT++)
for (size_t iD = 0; iD < dofs_per_cell; iD++)
dof_offsets[iT*dofs_per_cell+iD] = e.cell_global_dof_offset(iT)+iD;
}
};
inline std::ostream& operator<<(std::ostream& os, const dof_range& dr)
{
os << "Tag " << dr.tag << ", dim " << dr.dim << ": (";
os << dr.base << ", " << dr.dof_offsets.size() << ")";
return os;
}
#endif
class model
{
/* The objects of the class 'model' are always local to the processor
* and they contain only the entities pertaining to the specific
* processor. In the 'model' class however there are additional
* data structures, which are needed by Rank 0 to keep track of
* the whole GMSH unpartitioned model. In other words, the 'model'
* class reflects more a GMSH partition. This mismatch in terminology
* is due to the fact that MPI was added in a second time.
* In the code, the term 'model' will thus refer to the GMSH partition,
* and will be a single computational unit.
* The collection of the models will be referred as 'world'.
*/
int geometric_order;
int approximation_order;
#ifdef USE_MPI
int num_partitions;
#endif /* USE_MPI */
......@@ -141,10 +122,10 @@ class model
/* Map from 3D entity tag to partition */
std::map<int, int> partition_inverse_map;
public:
/* RANK 0 only: keeps global (process-wide) DOF numbering and positioning */
std::map<int, dof_range> global_dof_ranges;
private:
/* RANK 0 only: keeps track of all the existing entities */
std::vector<int> all_entities_tags;
/* RANK 0 only: copy of remote entities for postpro */
std::map<int, std::vector<entity>> remote_entities;
#endif /* USE_MPI */
......@@ -184,20 +165,25 @@ public:
void populate_from_gmsh(void);
#ifdef USE_MPI
const std::map<int, dof_range>& dof_ranges(void) const
const std::vector<int>& world_entities_tags(void) const
{
return global_dof_ranges;
return all_entities_tags;
}
size_t global_num_dofs() const
size_t num_dofs_world() const
{
size_t ret = 0;
for (auto& dr : global_dof_ranges)
ret += dr.second.dof_offsets.size();
for (auto& e : entities)
ret += e.num_dofs();
for (auto& re : remote_entities)
for (auto& e : re.second)
ret += e.num_dofs();
return ret;
}
size_t global_num_cells() const
size_t num_cells_world() const
{
size_t ret = 0;
for (auto& e : entities)
......@@ -215,8 +201,8 @@ public:
return partition_inverse_map.at(tag)-1;
}
entity& entity_global_lookup(int);
const entity& entity_global_lookup(int) const;
entity& entity_world_lookup(int);
const entity& entity_world_lookup(int) const;
#endif /* USE_MPI */
const element_connectivity<element_key>& connectivity() const
......
src/entity.cpp
+ 65
− 29
View file @ cdcc53a1
......@@ -9,7 +9,9 @@
entity::entity(const entity_params& ep)
: dim(ep.dim), tag(ep.tag), elemType(ep.etype), g_order(ep.gorder),
a_order(ep.aorder), m_dof_base(0), m_flux_base(0), m_index_base(0),
a_order(ep.aorder), m_dof_base_model(0), m_flux_base_model(0),
m_index_base_model(0), m_entity_number_model(0), m_dof_base_world(0),
m_flux_base_world(0), m_index_base_world(0), m_entity_number_world(0),
cur_elem_ordering(entity_ordering::GMSH)
{
/* Prepare reference elements */
......@@ -98,7 +100,7 @@ entity::project(const scalar_function& function) const
}
/* Project a function on this entity and store the partial result in
* the global vector pf at the offset corresponding to this entity. */
* the model-local vector pf at the offset corresponding to this entity. */
void
entity::project(const scalar_function& function, vecxd& pf) const
{
......@@ -124,12 +126,12 @@ entity::project(const scalar_function& function, vecxd& pf) const
rhs += pqp.weight() * function(pqp.point()) * phi;
}
pf.segment(m_dof_base + iT*num_bf, num_bf) = mass.ldlt().solve(rhs);
pf.segment(m_dof_base_model + iT*num_bf, num_bf) = mass.ldlt().solve(rhs);
}
}
/* Project a function on this entity and store the partial result in
* the global vector pf at the offset corresponding to this entity. */
* the model-local vector pf at the offset corresponding to this entity. */
void
entity::project(const vector_function& function, vecxd& pfx, vecxd& pfy, vecxd& pfz) const
{
......@@ -162,9 +164,9 @@ entity::project(const vector_function& function, vecxd& pfx, vecxd& pfy, vecxd&
Eigen::LDLT<matxd> mass_ldlt;
mass_ldlt.compute(mass);
pfx.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_x);
pfy.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_y);
pfz.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_z);
pfx.segment(m_dof_base_model + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_x);
pfy.segment(m_dof_base_model + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_y);
pfz.segment(m_dof_base_model + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_z);
}
}
......@@ -206,9 +208,9 @@ entity::project(const vector_function& function, double *pfx, double *pfy, doubl
for (size_t i = 0; i < num_bf; i++)
{
pfx[m_dof_base + iT*num_bf + i] = p_x[i];
pfy[m_dof_base + iT*num_bf + i] = p_y[i];
pfz[m_dof_base + iT*num_bf + i] = p_z[i];
pfx[m_dof_base_model + iT*num_bf + i] = p_x[i];
pfy[m_dof_base_model + iT*num_bf + i] = p_y[i];
pfz[m_dof_base_model + iT*num_bf + i] = p_z[i];
}
}
}
......@@ -267,10 +269,17 @@ entity::cell_refelem(const physical_element& pe) const
}
size_t
entity::cell_global_index(size_t iT) const
entity::cell_model_index(size_t iT) const
{
assert(iT < physical_cells.size());
return m_index_base + iT;
return m_index_base_model + iT;
}
size_t
entity::cell_world_index(size_t iT) const
{
assert(iT < physical_cells.size());
return m_index_base_world + iT;
}
size_t
......@@ -281,12 +290,18 @@ entity::cell_local_index_by_gmsh(size_t iT) const
}
size_t
entity::cell_global_index_by_gmsh(size_t iT) const
entity::cell_model_index_by_gmsh(size_t iT) const
{
assert(iT < physical_cells.size());
return m_index_base + physical_cells[iT].original_position();
return m_index_base_model + physical_cells[iT].original_position();
}
size_t
entity::cell_world_index_by_gmsh(size_t iT) const
{
assert(iT < physical_cells.size());
return m_index_base_world + physical_cells[iT].original_position();
}
size_t
entity::cell_local_dof_offset(size_t iT) const
......@@ -296,10 +311,17 @@ entity::cell_local_dof_offset(size_t iT) const
}
size_t
entity::cell_global_dof_offset(size_t iT) const
entity::cell_model_dof_offset(size_t iT) const
{
assert(reference_cells.size() > 0);
return m_dof_base_model + iT * reference_cells[0].num_basis_functions();
}
size_t
entity::cell_world_dof_offset(size_t iT) const
{
assert(reference_cells.size() > 0);
return m_dof_base + iT * reference_cells[0].num_basis_functions();
return m_dof_base_world + iT * reference_cells[0].num_basis_functions();
}
const reference_element&
......@@ -369,41 +391,55 @@ entity::num_fluxes(void) const
}
void
entity::base(size_t d_base, size_t f_base, size_t i_base)
entity::base_model(size_t d_base, size_t f_base, size_t i_base)
{
m_dof_base = d_base;
m_flux_base = f_base;
m_index_base = i_base;
m_dof_base_model = d_base;
m_flux_base_model = f_base;
m_index_base_model = i_base;
}
void
entity::base_world(size_t d_base, size_t f_base, size_t i_base)
{
m_dof_base_world = d_base;
m_flux_base_world = f_base;
m_index_base_world = i_base;
}
size_t
entity::dof_base(void) const
{
return m_dof_base;
return m_dof_base_model;
}
size_t
entity::dof_base_world(void) const
{
return m_dof_base_world;
}
size_t
entity::flux_base(void) const
{
return m_flux_base;
return m_flux_base_model;
}
size_t
entity::index_base(void) const
{
return m_index_base;
return m_index_base_model;
}
size_t
entity::number(void) const
{
return m_entity_number;
return m_entity_number_model;
}
void
entity::number(size_t n)
{
m_entity_number = n;
m_entity_number_model = n;
}
matxd
......@@ -848,8 +884,8 @@ entity::populate_entity_data(entity_data_cpu& ed, const model& mod) const
ed.num_fluxes = reference_faces[0].num_basis_functions();
ed.num_qp = (g_order == 1) ? 1 : reference_cells[0].num_integration_points();
ed.num_face_qp = (g_order == 1) ? 1 : reference_faces[0].num_integration_points();
ed.dof_base = m_dof_base;
ed.flux_base = m_flux_base;
ed.dof_base = m_dof_base_model;
ed.flux_base = m_flux_base_model;
ed.num_faces_per_elem = num_faces_per_elem();
populate_differentiation_matrices(ed);
......@@ -885,7 +921,7 @@ entity::populate_entity_data(entity_data_cpu& ed, const model& mod) const
for (size_t i = 0; i < face_keys.size(); i++)
for (size_t j = 0; j < cell_keys.size(); j++)
if ( face_keys.at(i) == cell_keys.at(j) )
ed.fluxdofs_mine.at(face_base + i) = cell_global_dof_offset(iT) + j;
ed.fluxdofs_mine.at(face_base + i) = cell_model_dof_offset(iT) + j;
element_key fk(pf);
......@@ -897,7 +933,7 @@ entity::populate_entity_data(entity_data_cpu& ed, const model& mod) const
for (size_t i = 0; i < face_keys.size(); i++)
for (size_t j = 0; j < ne_cell_keys.size(); j++)
if ( face_keys.at(i) == ne_cell_keys.at(j) )
ed.fluxdofs_neigh.at(face_base + i) = ne_e.cell_global_dof_offset(neigh.iT) + j;
ed.fluxdofs_neigh.at(face_base + i) = ne_e.cell_model_dof_offset(neigh.iT) + j;
}
}
}
......
src/gmsh_io.cpp
+ 52
− 38
View file @ cdcc53a1
......@@ -195,13 +195,14 @@ model::import_gmsh_entities_rank0(void)
gvp_t gmsh_entities;
gm::getEntities(gmsh_entities, DIMENSION(3));
size_t entity_index = 0;
size_t dof_base = 0;
size_t flux_base = 0;
size_t index_base = 0;
#ifdef USE_MPI
size_t global_dof_base = 0;
#endif /* USE_MPI */
size_t entity_index_model = 0;
size_t dof_base_model = 0;
size_t flux_base_model = 0;
size_t index_base_model = 0;
size_t entity_index_world = 0;
size_t dof_base_world = 0;
size_t flux_base_world = 0;
size_t index_base_world = 0;
for (auto [dim, tag] : gmsh_entities)
{
......@@ -245,16 +246,22 @@ model::import_gmsh_entities_rank0(void)
if (entity_is_remote)
{
dof_range dr(e, global_dof_base);
global_dof_ranges.emplace(tag, std::move(dr));
global_dof_base += e.num_dofs();
/* Update world-offsets */
e.base_world(dof_base_world, flux_base_world, index_base_world);
//e.number(entity_index_world);
dof_base_world += e.num_dofs();
flux_base_world += e.num_fluxes();
index_base_world += e.num_cells();
/* Send entity to the appropriate process */
assert(tag_partition > 0);
int rank = tag_partition - 1;
e.mpi_send(rank, MPI_COMM_WORLD);
std::cout << "Sending entity " << tag << " to " << rank << " " << e.num_dofs() << std::endl;
/* Save locally the entity */
remote_entities[tag_partition].push_back( std::move(e) );
all_entities_tags.push_back(tag);
entity_index_world++;
continue;
}
}
......@@ -267,24 +274,30 @@ model::import_gmsh_entities_rank0(void)
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);
etag_to_entity_offset[ pe.element_tag() ] = std::make_pair(entity_index_model, i);
}
e.base(dof_base, flux_base, index_base);
e.number(entity_index);
dof_base += e.num_dofs();
flux_base += e.num_fluxes();
index_base += e.num_cells();
/* Update model-offsets */
e.base_model(dof_base_model, flux_base_model, index_base_model);
e.number(entity_index_model);
dof_base_model += e.num_dofs();
flux_base_model += e.num_fluxes();
index_base_model += e.num_cells();
#ifdef USE_MPI
dof_range dr(e, global_dof_base);
global_dof_ranges.emplace(tag, std::move(dr));
global_dof_base += e.num_dofs();
#endif /* USE_MPI */
/* Update world-offsets */
e.base_world(dof_base_world, flux_base_world, index_base_world);
//e.number(entity_index_world);
dof_base_world += e.num_dofs();
flux_base_world += e.num_fluxes();
index_base_world += e.num_cells();
update_connectivity(e, entity_index);
update_connectivity(e, entity_index_model);
entities.push_back( std::move(e) );
entity_index++;
#ifdef USE_MPI
all_entities_tags.push_back(tag);
#endif /* USE_MPI */
entity_index_model++;
entity_index_world++;
}
}
......@@ -293,20 +306,21 @@ model::import_gmsh_entities_rank0(void)
assert( IMPLIES(num_partitions == 1, remote_entities.size() == 0) );
for (auto& re : remote_entities)
{
entity_index = 0;
dof_base = 0;
flux_base = 0;
index_base = 0;
entity_index_model = 0;
dof_base_model = 0;
flux_base_model = 0;
index_base_model = 0;
for (auto& e : re.second)
{
e.base(dof_base, flux_base, index_base);
e.number(entity_index);
dof_base += e.num_dofs();
flux_base += e.num_fluxes();
index_base += e.num_cells();
entity_index++;
e.base_model(dof_base_model, flux_base_model, index_base_model);
e.number(entity_index_model);
dof_base_model += e.num_dofs();
flux_base_model += e.num_fluxes();
index_base_model += e.num_cells();
entity_index_model++;
}
}
std::sort(all_entities_tags.begin(), all_entities_tags.end());
#endif
}
......@@ -337,7 +351,7 @@ model::import_gmsh_entities_rankN(void)
etag_to_entity_offset[ pe.element_tag() ] = std::make_pair(entity_index, i);
}
e.base(dof_base, flux_base, index_base);
e.base_model(dof_base, flux_base, index_base);
e.number(entity_index);
dof_base += e.num_dofs();
flux_base += e.num_fluxes();
......@@ -585,7 +599,7 @@ model::partition_mesh(int parts)
}
entity&
model::entity_global_lookup(int tag)
model::entity_world_lookup(int tag)
{
/* Look for tag between local entities */
for (auto& e : entities)
......@@ -602,7 +616,7 @@ model::entity_global_lookup(int tag)
}
const entity&
model::entity_global_lookup(int tag) const
model::entity_world_lookup(int tag) const
{
/* Look for tag between local entities */
for (const auto& e : entities)
......
src/maxwell_common.cpp
+ 9
− 9
View file @ cdcc53a1
......@@ -156,7 +156,7 @@ void export_E_field_nodal(const model& mod, silo& sdb, const vecxd& Ex,
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
auto cgofs = e.cell_global_dof_offset(iT);
auto cgofs = e.cell_model_dof_offset(iT);
auto nodetags = pe.node_tags();
for (size_t iD = 0; iD < 4; iD++)
{
......@@ -196,7 +196,7 @@ void export_H_field_nodal(const model& mod, silo& sdb, const vecxd& Hx,
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
auto cgofs = e.cell_global_dof_offset(iT);
auto cgofs = e.cell_model_dof_offset(iT);
auto nodetags = pe.node_tags();
for (size_t iD = 0; iD < 4; iD++)
{
......@@ -239,7 +239,7 @@ void export_J_field_nodal(const model& mod, silo& sdb, const vecxd& Ex,
auto& pe = e.cell(iT);
auto bar = pe.barycenter();
auto sigma = mpl.sigma(etag, bar);
auto cgofs = e.cell_global_dof_offset(iT);
auto cgofs = e.cell_model_dof_offset(iT);
auto nodetags = pe.node_tags();
for (size_t iD = 0; iD < 4; iD++)
{
......@@ -284,8 +284,8 @@ void export_E_field_zonal(const model& mod, silo& sdb, const vecxd& Ex,
auto& re = e.cell_refelem(pe);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
auto num_bf = re.num_basis_functions();
auto ofs = e.cell_global_dof_offset(iT);
auto gi = e.cell_global_index_by_gmsh(iT);
auto ofs = e.cell_model_dof_offset(iT);
auto gi = e.cell_model_index_by_gmsh(iT);
out_Ex[gi] = Ex.segment(ofs, num_bf).dot(phi_bar);
out_Ey[gi] = Ey.segment(ofs, num_bf).dot(phi_bar);
out_Ez[gi] = Ez.segment(ofs, num_bf).dot(phi_bar);
......@@ -314,8 +314,8 @@ void export_H_field_zonal(const model& mod, silo& sdb, const vecxd& Hx,
auto& re = e.cell_refelem(pe);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
auto num_bf = re.num_basis_functions();
auto ofs = e.cell_global_dof_offset(iT);
auto gi = e.cell_global_index_by_gmsh(iT);
auto ofs = e.cell_model_dof_offset(iT);
auto gi = e.cell_model_index_by_gmsh(iT);
out_Hx[gi] = Hx.segment(ofs, num_bf).dot(phi_bar);
out_Hy[gi] = Hy.segment(ofs, num_bf).dot(phi_bar);
out_Hz[gi] = Hz.segment(ofs, num_bf).dot(phi_bar);
......@@ -347,8 +347,8 @@ void export_J_field_zonal(const model& mod, silo& sdb, const vecxd& Ex,
auto sigma = mpl.sigma(etag, bar);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
auto num_bf = re.num_basis_functions();
auto ofs = e.cell_global_dof_offset(iT);
auto gi = e.cell_global_index_by_gmsh(iT);
auto ofs = e.cell_model_dof_offset(iT);
auto gi = e.cell_model_index_by_gmsh(iT);
out_Jx[gi] = sigma*Ex.segment(ofs, num_bf).dot(phi_bar);
out_Jy[gi] = sigma*Ey.segment(ofs, num_bf).dot(phi_bar);
out_Jz[gi] = sigma*Ez.segment(ofs, num_bf).dot(phi_bar);
......
src/maxwell_cpu.cpp
+ 24
− 30
View file @ cdcc53a1
......@@ -31,7 +31,7 @@ init_matparams(const model& mod, solver_state& state,
auto epsilon = mpl.epsilon(etag, bar);
auto mu = mpl.mu(etag, bar);
auto sigma = mpl.sigma(etag, bar);
auto ofs = e.cell_global_dof_offset(iT);
auto ofs = e.cell_model_dof_offset(iT);
for (size_t iD = 0; iD < re.num_basis_functions(); iD++)
{
......@@ -600,39 +600,37 @@ gather_field(const model& mod, maxwell::field& in, maxwell::field& out,
if (rank == root)
{
auto& gdr = mod.dof_ranges();
for (auto& dr : gdr)
auto& all_tags = mod.world_entities_tags();
for (auto& tag : all_tags)
{
auto tag = dr.second.tag;
if (mod.entity_rank(tag) == root)
continue;
auto base = dr.second.base;
auto len = dr.second.dof_offsets.size();
//std::cout << "GR: " << tag << " " << len << std::endl;
auto& e = mod.entity_world_lookup(tag);
auto base = e.dof_base_world();
auto len = e.num_dofs();
std::cout << "GR: " << tag << " " << len << std::endl;
receive_field(out, mod.entity_rank(tag), tag, base, len, comm);
}
for (auto& e : mod)
{
auto in_base = e.dof_base();
auto in_len = e.num_dofs();
auto out_base = gdr.at(e.gmsh_tag()).base;
auto out_len = gdr.at(e.gmsh_tag()).dof_offsets.size();
assert(in_len == out_len);
out.Ex.segment(out_base, out_len) = in.Ex.segment(in_base, in_len);
out.Ey.segment(out_base, out_len) = in.Ey.segment(in_base, in_len);
out.Ez.segment(out_base, out_len) = in.Ez.segment(in_base, in_len);
out.Hx.segment(out_base, out_len) = in.Hx.segment(in_base, in_len);
out.Hy.segment(out_base, out_len) = in.Hy.segment(in_base, in_len);
out.Hz.segment(out_base, out_len) = in.Hz.segment(in_base, in_len);
auto len = e.num_dofs();
auto out_base = e.dof_base_world();
out.Ex.segment(out_base, len) = in.Ex.segment(in_base, len);
out.Ey.segment(out_base, len) = in.Ey.segment(in_base, len);
out.Ez.segment(out_base, len) = in.Ez.segment(in_base, len);
out.Hx.segment(out_base, len) = in.Hx.segment(in_base, len);
out.Hy.segment(out_base, len) = in.Hy.segment(in_base, len);
out.Hz.segment(out_base, len) = in.Hz.segment(in_base, len);
}
}
else
{
for (auto& e : mod)
{
//std::cout << "GS: " << e.gmsh_tag() << " " << e.num_dofs() << std::endl;
std::cout << "GS: " << e.gmsh_tag() << " " << e.num_dofs() << std::endl;
send_field(in, root, e.gmsh_tag(), e.dof_base(), e.num_dofs(), comm);
}
}
......@@ -648,7 +646,7 @@ export_fields_to_silo_2(const model& mod, maxwell::solver_state& state,
if (rank == 0)
{
size_t gdnum = mod.global_num_dofs();
size_t gdnum = mod.num_dofs_world();
maxwell::field f;
f.resize(gdnum);
gather_field(mod, state.emf_next, f, 0, MPI_COMM_WORLD);
......@@ -661,29 +659,25 @@ export_fields_to_silo_2(const model& mod, maxwell::solver_state& state,
sdb.import_mesh_from_gmsh();
sdb.write_mesh(state.curr_time, state.curr_timestep);
std::vector<double> out_Ex( mod.global_num_cells() );
std::vector<double> out_Ey( mod.global_num_cells() );
std::vector<double> out_Ez( mod.global_num_cells() );
std::vector<double> out_Ex( mod.num_cells_world() );
std::vector<double> out_Ey( mod.num_cells_world() );
std::vector<double> out_Ez( mod.num_cells_world() );
size_t nc = 0;
for (auto& dr : mod.dof_ranges())
for (auto& tag : mod.world_entities_tags())
{
auto etag = dr.first;
std::cout << etag << std::endl;
auto& e = mod.entity_global_lookup(etag);
auto& e = mod.entity_world_lookup(tag);
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
auto& re = e.cell_refelem(pe);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
auto num_bf = re.num_basis_functions();
auto ofs = dr.second.base + num_bf*iT;
auto gi = nc+e.cell_local_index_by_gmsh(iT);
auto ofs = e.cell_world_dof_offset(iT);
auto gi = e.cell_world_index_by_gmsh(iT);
out_Ex[gi] = f.Ex.segment(ofs, num_bf).dot(phi_bar);
out_Ey[gi] = f.Ey.segment(ofs, num_bf).dot(phi_bar);
out_Ez[gi] = f.Ez.segment(ofs, num_bf).dot(phi_bar);
}
nc += e.num_cells();
}
sdb.write_zonal_variable("Ex", out_Ex);
......
src/maxwell_solver.cpp
+ 2
− 2
View file @ cdcc53a1
......@@ -139,7 +139,7 @@ integrate_electric_field(const model& mod, const maxwell::solver_state& state,
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_ofs = e.cell_model_dof_offset(ofs);
auto dof_num = re.num_basis_functions();
vecxd phi = re.basis_functions( refp );
vecxd locEx = state.emf_curr.Ex.segment(dof_ofs, dof_num);
......@@ -175,7 +175,7 @@ reinit_materials(const model& mod, maxwell::solver_state& state, maxwell::parame
auto epsilon = mpl.epsilon(tag, bar);
//auto mu = mpl.mu(tag, bar);
auto sigma = mpl.sigma(tag, bar);
auto ofs = e.cell_global_dof_offset(iT);
auto ofs = e.cell_model_dof_offset(iT);
for (size_t iD = 0; iD < re.num_basis_functions(); iD++)
{
......
tests/test_differentiation.cpp
+ 1
− 1
View file @ cdcc53a1
......@@ -147,7 +147,7 @@ int test_differentiation_convergence(int geometric_order, int approximation_orde
auto& re = e.cell_refelem(pe);
matxd mass = e.mass_matrix(iT);
auto num_bf = re.num_basis_functions();
auto ofs = e.cell_global_dof_offset(iT);
auto ofs = e.cell_model_dof_offset(iT);
vecxd diff_x = Pdf_dx.segment(ofs, num_bf) - Cdf_dx.segment(ofs, num_bf);
vecxd diff_y = Pdf_dy.segment(ofs, num_bf) - Cdf_dy.segment(ofs, num_bf);
vecxd diff_z = Pdf_dz.segment(ofs, num_bf) - Cdf_dz.segment(ofs, num_bf);
......
tests/test_differentiation_gpu.cpp
+ 1
− 1
View file @ cdcc53a1
......@@ -205,7 +205,7 @@ int test_differentiation_convergence(int geometric_order, int approximation_orde
auto& re = e.cell_refelem(pe);
matxd mass = e.mass_matrix(iT);
auto num_bf = re.num_basis_functions();
auto ofs = e.cell_global_dof_offset(iT);
auto ofs = e.cell_model_dof_offset(iT);
vecxd diff_x = Pdf_dx.segment(ofs, num_bf) - Cdf_dx.segment(ofs, num_bf);
vecxd diff_y = Pdf_dy.segment(ofs, num_bf) - Cdf_dy.segment(ofs, num_bf);
vecxd diff_z = Pdf_dz.segment(ofs, num_bf) - Cdf_dz.segment(ofs, num_bf);
......
tests/test_lifting.cpp
+ 1
− 1
View file @ cdcc53a1
......@@ -198,7 +198,7 @@ int test_lifting(int geometric_order, int approximation_order)
vol += pqp.weight() * (dphi*J.inverse()) * F(pqp.point());
}
auto ofs = e.cell_global_dof_offset(iT);
auto ofs = e.cell_model_dof_offset(iT);
vecxd excpected_div_F = Pdiv_F.segment(ofs, num_bf);
vecxd Plf = LiftF.segment(ofs, num_bf);
vecxd Pvol = mass.ldlt().solve(vol);
......
tests/test_lifting_gpu.cpp
+ 1
− 1
View file @ cdcc53a1
......@@ -207,7 +207,7 @@ int test_lifting(int geometric_order, int approximation_order)
vol += pqp.weight() * (dphi*J.inverse()) * F(pqp.point());
}
auto ofs = e.cell_global_dof_offset(iT);
auto ofs = e.cell_model_dof_offset(iT);
vecxd excpected_div_F = Pdiv_F.segment(ofs, num_bf);
vecxd Plf = LiftF.segment(ofs, num_bf);
vecxd Pvol = mass.ldlt().solve(vol);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment