Merge pull request #514 from sblauth/dev/mpi_topology

Fix bugs for topology optimization in parallel
sblauth · Nov 11, 2024 · cba7928 · cba7928
2 parents a434638 + 902f37a
commit cba7928
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 3 deletions.
diff --git a/cashocs/_optimization/topology_optimization/topology_optimization_algorithm.py b/cashocs/_optimization/topology_optimization/topology_optimization_algorithm.py
@@ -90,6 +90,16 @@ def __init__(
         self._cashocs_problem.db.parameter_db.problem_type = "topology"
 
         self.mesh = optimization_problem.mesh
+
+        mpi_comm = self.mesh.mpi_comm()
+        if self.interpolation_scheme == "angle" and mpi_comm.Get_size() > 1:
+            raise _exceptions.InputError(
+                "TopologyOptimizationProblem",
+                "TopologyOptimization.interpolation_scheme",
+                "The angle weighted interpolation option is not supported in parallel. "
+                "Please use interpolation_scheme = volume in your config file.",
+            )
+
         self.cg1_space = fenics.FunctionSpace(self.mesh, "CG", 1)
         self.dg0_space = optimization_problem.dg0_space
         self.topological_derivative_vertex: fenics.Function = fenics.Function(

diff --git a/cashocs/_utils/interpolations.py b/cashocs/_utils/interpolations.py
@@ -132,10 +132,11 @@
 
     std::vector<unsigned int> cells = mesh->cells();
     int meshdim = mesh->geometry().dim();
+    auto ghost_offset = mesh->topology().ghost_offset(meshdim);
 
     if (meshdim == 2){
         int index = 0;
-        for (int i=0; i<cells.size(); i+=3){
+        for (int i=0; i<ghost_offset*3; i+=3){
             psi0 = vertex_values[cells[i]];
             psi1 = vertex_values[cells[i+1]];
             psi2 = vertex_values[cells[i+2]];
@@ -155,7 +156,7 @@
     }
     else if (meshdim == 3){
         int index = 0;
-        for (int i=0; i<cells.size(); i+=4){
+        for (int i=0; i<ghost_offset*4; i+=4){
             psi0 = vertex_values[cells[i]];
             psi1 = vertex_values[cells[i+1]];
             psi2 = vertex_values[cells[i+2]];
@@ -246,7 +247,8 @@ def interpolate_by_volume(
     )
     arr = fenics.assemble(form_td * test * dx)
     vol = fenics.assemble(test * dx)
-    node_function.vector()[:] = arr[:] / vol[:]
+    node_function.vector().set_local(arr[:] / vol[:])
+    node_function.vector().apply("")
 
 
 def interpolate_by_angle(
@@ -439,3 +441,4 @@ def interpolate_by_angle(
     values /= weights
     d2v = fenics.dof_to_vertex_map(cg1_space)
     node_function.vector()[:] = values[d2v]
+    node_function.vector().apply("")