clean _init_dof_maps

LDDsimulation/domainPatch.py

@@ -379,9 +379,9 @@ class DomainPatch(df.SubDomain):
                         # assumes a Richards model, we assume constant normalised atmospheric pressure
                         # and no interface term is practically appearing.
                         # interface_forms += (df.Constant(0)*phi_a)*ds(interface)
-                        # pass
-                        interface_forms.append((0*pa*phi_a)*ds(marker_value))
-                        gli_forms.append((0*phi_a)*ds(marker_value))
+                        pass
+                        # interface_forms.append((0*pa*phi_a)*ds(marker_value))
+                        # gli_forms.append((0*phi_a)*ds(marker_value))
                     else:
                         # interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
                         interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
@@ -797,47 +797,35 @@ class DomainPatch(df.SubDomain):
             # we want the same space for both primary variables (pw
             # and pnw) to be able to calculate pnw-pw without issues.
             pressure_space = df.FunctionSpace(self.mesh, 'P', degree)
+            gli_space = df.FunctionSpace(self.mesh, 'DG', degree)
+            flux_space = df.VectorFunctionSpace(self.mesh, 'DG', degree)
             for phase in self.has_phases:
                 if function == "pressure":
                     self.function_space[function].update({phase: pressure_space})
                     self.trialfunction.update({phase: df.TrialFunction(self.function_space["pressure"][phase])})
                     self.testfunction.update({phase: df.TestFunction(self.function_space["pressure"][phase])})
                 elif function == "gli":
-                    degree = self.function_space["pressure"][phase].ufl_element().degree()
-                    self.function_space[function].update({phase: df.FunctionSpace(self.mesh, 'DG', degree)})
+                    self.function_space[function].update({phase: gli_space})
                 else:
                     # here function = flux
-                    degree = self.function_space["pressure"][phase].ufl_element().degree()
-                    self.function_space[function].update({phase: df.VectorFunctionSpace(self.mesh, 'DG', degree)})
+                    # degree = self.function_space["pressure"][phase].ufl_element().degree()
+                    self.function_space[function].update({phase: flux_space})
 
 
     def _init_dof_maps(self) -> None:
-        """ calculate dof maps and the vertex to parent map.
+        """ calculate and save dof maps of the different function spaces. This is needed for
+        writing the communication over the interface.
 
-        This method sets the class Variables
-            self.dof2vertex
-            self.vertex2dof
+        This method populates the class dictionaries
+            self.dofmap
+        Note, that the flux dofmap is split into x and y coordinates.
         """
         mesh_data = self.mesh.data()
-        # # local to parent vertex index map
-        # if self.has_interface is not None:
-        #     self.parent_mesh_index = mesh_data.array('parent_vertex_indices',0)
-        # print(f"on subdomain{self.subdomain_index} parent_vertex_indices: {self.parent_mesh_index}")
-        self.dof2vertex = dict()
-        self.vertex2dof = dict()
         self.dofmap = dict()
         self.dofmap.update({"pressure": dict()})
         self.dofmap.update({"gli": dict()})
         self.dofmap.update({"flux": dict()})
         for phase in self.has_phases:
-            self.dof2vertex.update(#
-                {phase :df.dof_to_vertex_map(self.function_space["pressure"][phase])}#
-                )
-            # print(f"on subdomain{self.subdomain_index} dof2vertex: {self.dof2vertex[phase]}")
-            self.vertex2dof.update(#
-                {phase :df.vertex_to_dof_map(self.function_space["pressure"][phase])}#
-                )
-            # print(f"on subdomain{self.subdomain_index} vertex2dof: {self.vertex2dof[phase]}")
             self.dofmap["pressure"].update(
                 {phase: self.function_space["pressure"][phase].dofmap()}
                 )
@@ -1125,9 +1113,12 @@ class DomainPatch(df.SubDomain):
         for phase in self.has_phases:
             rho = df.Constant(self.densities[phase])
             self.gravity_term.update(
-                {phase: df.Expression('-g*rho*x[1]', domain = self.mesh, #
+                {phase: df.Expression('-g*rho*x[1]',
+                                      domain=self.mesh, #
                                       degree=self.interpolation_degree,
-                                       g = g, rho = rho)}
+                                      g=g,
+                                      rho=rho
+                                      )}
             )
 
     def write_solution_to_xdmf(self, file: tp.Type[SolutionFile], #