diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py index f2bbad117cc1ff67bef82c0e01671489634b98ca..92d2aef7de3fe16f2dcb528aba2f9486f871fc6a 100644 --- a/LDDsimulation/domainPatch.py +++ b/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, # - degree = self.interpolation_degree, - g = g, rho = rho)} + {phase: df.Expression('-g*rho*x[1]', + domain=self.mesh, # + degree=self.interpolation_degree, + g=g, + rho=rho + )} ) def write_solution_to_xdmf(self, file: tp.Type[SolutionFile], #