Skip to content
Snippets Groups Projects
Commit b4aa88df authored by David Seus's avatar David Seus
Browse files

fix wrong assignment of gli_dofs on interface facets

parent ed4e464b
No related branches found
No related tags found
No related merge requests found
......@@ -57,6 +57,8 @@
*.h5
*.dep
*.csv
*.ipynb
core
# Ignoriere Bilder und Graphiken sowie Videos und Musik
*.png
......
This diff is collapsed.
......@@ -333,7 +333,8 @@ class Interface(BoundaryPart):
interface_marker_value: int,
debug: bool = False) -> tp.Dict[int, tp.Dict[int, np.ndarray]]:
""" for a given function space function_space, calculate the dof indices
and coordinates along the interface.
and coordinates along the interface. The pair {global_dof_index: dof_coordinate}
is saved as value in a dictionary which has global facet indices as keys.
# Parameters
function_space: #type df.functionSpace the function space of whitch
......@@ -362,9 +363,21 @@ class Interface(BoundaryPart):
for mesh_entity in mesh_entity_iterator:
# this is the facet index relative to the Mesh
# that interface_marker is defined on.
facet_index = mesh_entity.index()
interface_facet_index = mesh_entity.index()
for interface_facet in df.facets(mesh_entity):
# we only need the x,y coordinates. vertex.point()[:] returns
# the array and cannot be sliced. the result can be indexed, hence,
# vertex.point()[:][0:2]
# facet_midpoint = interface_facet.midpoint()[:][0:2]
# print(f"facet: {interface_facet}")
# print("facet coordinates: ")
facet_points = []
for vertex in df.vertices(interface_facet):
# print(f"vertex: {vertex} with coordinates {vertex.point()[:][0:2]}")
facet_points.append(vertex.point())
dofs_and_coords.update({facet_index: dict()})
dofs_and_coords.update({interface_facet_index: dict()})
# because mesh_entity doesn't have a method to access coordinates we
# need to cast it into an acutal facet object.
for cell in df.cells(mesh_entity):
......@@ -373,12 +386,15 @@ class Interface(BoundaryPart):
# returns an iterator object.
facet_cell = cell
facet_cell_index = facet_cell.index()
# cell_facets_entity = df.facets(cell)
# cell_coordinates = cell.get_vertex_coordinates()
# the cell is a triangle and for each dof numbered locally for that
# cell, we have the coordinates in a numpy array.
# cell, we have the coordinates in a numpy array. The numbering is the
# numbering local to the cell.
facet_cell_dof_coordinates = element.tabulate_dof_coordinates(facet_cell)
# the cell dof map gives to each local (to the cell) numbering of the
# dofs the index of the dof realitve to the global dof numbering of
# dofs the index of the dofs realitve to the global dof numbering of
# space, i.e. the index to access the dof with function.vector()[index].
facet_cell_dofmap = dofmap.cell_dofs(facet_cell_index)
......@@ -387,9 +403,9 @@ class Interface(BoundaryPart):
# cell is a triangle and only has one facet that belongs to
# the interface. we want to store only dofs that lie on that particular
# facet
if self._is_on_boundary_part(dof_coord):
if self._is_on_line_segment(facet_points[0], facet_points[1], dof_coord):
global_dof_ind = facet_cell_dofmap[local_dof_ind]
dofs_and_coords[facet_index].update({global_dof_ind: dof_coord})
dofs_and_coords[interface_facet_index].update({global_dof_ind: dof_coord})
if debug:
......@@ -482,10 +498,11 @@ class Interface(BoundaryPart):
for local_facet_index, local_facet_vertex in subdom_facet2coord.items():
for global_facet_index, global_facet_vertex in global_facet2coord.items():
# vertices are not necessarily in the same order
vertices_are_equal = np.array_equal(local_facet_vertex[0], global_facet_vertex[0])
vertices_are_equal += np.array_equal(local_facet_vertex[0], global_facet_vertex[1])
vertices_are_equal += np.array_equal(local_facet_vertex[1], global_facet_vertex[0])
vertices_are_equal += np.array_equal(local_facet_vertex[1], global_facet_vertex[1])
# np.array_equal(local_facet_vertex[0], global_facet_vertex[0])
vertices_are_equal = np.allclose(local_facet_vertex[0], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
vertices_are_equal += np.allclose(local_facet_vertex[0], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
both_vertices_are_equal = (vertices_are_equal == 2)
if both_vertices_are_equal:
self.local_to_global_facet_map[subdomain_index].update(
......
This diff is collapsed.
......@@ -8,6 +8,9 @@ import sys
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import typing as tp
import sympy as sym
# from termcolor import colored
# from colormap import Colormap
#import fenicstools as ft
......@@ -18,174 +21,112 @@ def print_once(*args,**kwargs):
print(*args,**kwargs)
sys.stdout.flush()
# def fit_IC(mesh,u):
# """
# Fit given data u on mesh to certain functions.
# Useful in the case of a single droplet.
# Returns the fitted data as class for IC.
# """
#
# # get mesh coordiantes
# xx = mesh.coordinates()[:,0]
# yy = mesh.coordinates()[:,1]
# nx = np.size(xx)
# ny = np.size(yy)
#
# # get function values
# rho_data = np.zeros(nx)
# phase_field_data = np.zeros(nx)
# mu_data = np.zeros(nx)
#
# for i in range(0,nx):
# rho_data[i] = u.sub(0)(xx[i],yy[i])
# mu_data[i] = u.sub(4)(xx[i],yy[i])
# phase_field_data[i] = u.sub(2)(xx[i],yy[i])
#
# # define functions to fit
# def rho(x,a,b,cx,cy,dx,dy):
# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
# return val
#
# def phase_field(x,a,b,cx,cy,dx,dy):
# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
# return val
#
# def mu(x,a,b,cx,cy,dx,dy):
# val = a*np.tanh(b*np.sqrt((x[0]-cx)**2+(x[1]-cy)**2)+dx) + dy
# return val
#
# xx1 = np.concatenate((xx[:,np.newaxis],yy[:,np.newaxis]),axis= 1)
# # fit data to function
# popt_rho = curve_fit(rho, xx1.T , rho_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
# [a,b,cx,cy,dx,dy] = popt_rho[0]
# # print(a,b,cx,cy,dx,dy)
#
# popt_phase_field = curve_fit(phase_field, xx1.T , phase_field_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
# [a1,b1,cx1,cy1,dx1,dy1] = popt_phase_field[0]
# # print(a1,b1,cx1,cy1,dx1,dy1)
#
# popt_mu = curve_fit(mu, xx1.T , mu_data, p0=[1./np.pi,-10,.5,.3,1.,.5])
# [a2,b2,cx2,cy2,dx2,dy2] = popt_mu[0]
# # print(a2,b2,cx2,cy2,dx2,dy2)
#
# # plot result
# # res = 500;
# # [X,Y] = np.meshgrid(np.linspace(0.,1.,res),np.linspace(0.,1.,res))
# # Z = np.zeros((res,res))
# # for i in range(0,res):
# # for j in range(0,res):
# # Z[i,j] = phase_field([X[i,j],Y[i,j]],a1,b1,cx1,cy1,d1)
# # plt.imshow(Z,cmap = "viridis", origin ="lower")
# # plt.colorbar()
# # plt.show()
#
# class fitted_IC(Expression):
# def eval(self, values, x):
# # rho
# values[0] = rho(x,a,b,cx,cy,dx,dy)
#
# # velocity
# for i in range(2):
# values[i + 1] = 0.
# # c
# values[2 + 1] = phase_field(x,a1,b1,cx1,cy1,dx1,dy1)
# # tau
# values[2 + 2] = 0.0
# # mu
# values[2 + 3] = mu(x,a2,b2,cx2,cy2,dx2,dy2)
# # sigma
# for i in range(2):
# values[2 + 4 + i] = 0.
#
# def value_shape(self):
# return (4 + 2*2,)
#
# return fitted_IC
#
#
# def tensor_jump(v, n):
# """
# DG jump operator for vector valued functions
# """
# return df.outer(v("+"), n("+")) + df.outer(v("-"), n("-"))
#
# def print_constants(self,file=None):
# """
# Print the names and values of all df.Constant attributes.
# """
# for key in self.__dict__.keys():
# attribute = self.__dict__[key]
# if type(attribute) is type(df.Constant(0.)):
# print_once(attribute.name() + " = " + str(attribute.values()),file=file)
# for key in self.__dict__["_EOS"].__dict__:
# attribute = self.__dict__["_EOS"].__dict__[key]
# if type(attribute) is type(df.Constant(0.)):
# print_once(attribute.name() + " = " + str(attribute.values()),file=file)
#
# def get_interface_thickness(solution,ME,xmin=-.5,xmax=1.5,ymin=0.,ymax=1.,ny=1000):
# """
# Compute interface thickness from simulation result.
# The interface thickness is here defined as the region where 0.1 < c < 0.9.
# Parameter: solution, function space ME
# Meshrelated: xmin,xmax,ymin,ymax
# Number of generated samples: ny
# """
# xmid = .5*(xmin+xmax)
# yy = np.linspace(ymin,ymax,ny)
# xx = xmid*np.ones_like(yy)
# # evaluate phasefield along following points
# zz = np.column_stack((xx,yy))
# p = ft.Probes(zz.flatten(), ME.sub(2).collapse())
# p(solution.sub(2))
# # values collected on root process
# values= p.array(root=0)
# if df.MPI.rank(df.MPI.comm_world) == 0:
# try:
# # get interface_start
# interface_start = np.min(np.where(abs(values-0.1)<5.e-3))
# # get interface_end
# interface_end = np.min(np.where(abs(values-0.9)<5.e-3))
# # compute interface_thickness
# interface_thickness = abs(yy[interface_end]-yy[interface_start])
# return interface_thickness
# except:
# print_once(colored("Phasefield values out of range.","red"))
# return -1
#
#
# def plot_solution(u):
# """ Plot solution for poster or the like..."""
# # Define "CD"-Colormap
# red = [0.12549019607843137, 0.3843137254901961, 0.6431372549019608, 0.9019607843137255]
# blue = [0.6823529411764706, 0.788235294117647, 0.8941176470588236, 1.0]
# green = [0.33725490196078434, 0.5254901960784314, 0.7137254901960784, 0.9019607843137255]
# c = Colormap()
# mycmap = c.cmap( {'red':red[::-1], 'green':green[::-1], 'blue':blue[::-1]})
# matplotlib.cm.register_cmap(name='mycmap', cmap=mycmap)
# # plt.rcParams['image.cmap'] = 'mycmap'
# plt.rcParams['image.cmap'] = 'RdYlBu'
# matplotlib.rcParams.update({'font.size': 17})
# matplotlib.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
# matplotlib.rc('text', usetex=True)
#
# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
# # fixes to remove white contour lines
# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
# cnt = plot(u,mode="contourf",levels=256,antialiased = True,linestyles=None)
#
# for c in cnt.collections:
# c.set_edgecolor("face")
#
# plt.xlim(0.45, 0.8)
# plt.ylim(0.2, 0.55)
# plt.axis("off")
# plt.annotate("$+$ phase",color='w',xy=(0.5, .45))
# plt.annotate("$-$ phase",xy=(0.7, .3))
# plt.annotate("$\Omega$",xy=(0.75, .5),fontsize=25)
# plt.annotate("$\gamma$",xy=(0.647, .358),fontsize=25)
# plt.annotate("",xy=(0.62,0.375),xytext=(0.66,0.33),arrowprops=dict(arrowstyle="<->",
# connectionstyle="arc3",facecolor='black',linewidth=1))
# #plt.show()
# plt.tight_layout()
# plt.savefig('diffuse_simulation.pdf',bbox_inches='tight')
# return 0
def generate_exact_solution_expressions(
symbols: tp.Dict[int, int],
isRichards: tp.Dict[int, bool],
symbolic_pressure: tp.Dict[int, tp.Dict[str, str]],
symbolic_capillary_pressure: tp.Dict[int, str],
viscosity: tp.Dict[int, tp.List[float]],#
porosity: tp.Dict[int, tp.Dict[str, float]],#
relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
relative_permeability_prime: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],
saturation_pressure_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
saturation_pressure_relationship_prime: tp.Dict[int, tp.Callable[...,None]] = None,#
symbolic_S_pc_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
symbolic_S_pc_relationship_prime: tp.Dict[int, tp.Callable[...,None]] = None,#
densities: tp.Dict[int, tp.Dict[str, float]] = None,#
gravity_acceleration: float = 9.81,
include_gravity: bool = False,
) -> tp.Dict[str, tp.Dict[int, tp.Dict[str, str]] ]:
""" Generate exact solution, source and initial condition code from symbolic expressions"""
output = {
'source': dict(),
'initial_condition': dict(),
'exact_solution': dict()
}
x = symbols["x"]
y = symbols["y"]
t = symbols["t"]
# turn above symbolic code into exact solution for dolphin and
# construct the rhs that matches the above exact solution.
dtS = dict()
div_flux = dict()
if saturation_pressure_relationship is not None:
S_pc_sym = saturation_pressure_relationship
S_pc_sym_prime = saturation_pressure_relationship_prime
for subdomain, isR in isRichards.items():
dtS.update({subdomain: dict()})
div_flux.update({subdomain: dict()})
output['source'].update({subdomain: dict()})
output['exact_solution'].update({subdomain: dict()})
output['initial_condition'].update({subdomain: dict()})
if isR:
subdomain_has_phases = ["wetting"]
else:
subdomain_has_phases = ["wetting", "nonwetting"]
# conditional for S_pc_prime
pc = symbolic_capillary_pressure[subdomain]
dtpc = sym.diff(pc, t, 1)
dxpc = sym.diff(pc, x, 1)
dypc = sym.diff(pc, y, 1)
if saturation_pressure_relationship is not None:
S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
else:
S = symbolic_S_pc_relationship[subdomain]
dS = symbolic_S_pc_relationship_prime[subdomain]
for phase in subdomain_has_phases:
# Turn above symbolic expression for exact solution into c code
output['exact_solution'][subdomain].update(
{phase: sym.printing.ccode(symbolic_pressure[subdomain][phase])}
)
# save the c code for initial conditions
output['initial_condition'][subdomain].update(
{phase: sym.printing.ccode(symbolic_pressure[subdomain][phase].subs(t, 0))}
)
if phase == "nonwetting":
dtS[subdomain].update(
{phase: -porosity[subdomain]*dS*dtpc}
)
else:
dtS[subdomain].update(
{phase: porosity[subdomain]*dS*dtpc}
)
pa = symbolic_pressure[subdomain][phase]
dxpa = sym.diff(pa, x, 1)
dxdxpa = sym.diff(pa, x, 2)
dypa = sym.diff(pa, y, 1)
dydypa = sym.diff(pa, y, 2)
mu = viscosity[subdomain][phase]
ka = relative_permeability[subdomain][phase]
dka = relative_permeability_prime[subdomain][phase]
rho = densities[subdomain][phase]
g = gravity_acceleration
if phase == "nonwetting":
# x part of div(flux) for nonwetting
dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
# y part of div(flux) for nonwetting
if include_gravity:
dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(1-S)
else:
dydyflux = -1/mu*dka(1-S)*dS*dypc*dypa + 1/mu*dydypa*ka(1-S)
else:
# x part of div(flux) for wetting
dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
# y part of div(flux) for wetting
if include_gravity:
dydyflux = 1/mu*dka(S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(S)
else:
dydyflux = 1/mu*dka(S)*dS*dypc*(dypa) + 1/mu*dydypa*ka(S)
div_flux[subdomain].update({phase: dxdxflux + dydyflux})
contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
output['source'][subdomain].update(
{phase: sym.printing.ccode(contructed_rhs)}
)
return output
......@@ -9,5 +9,5 @@ class SolutionFile(df.XDMFFile):
df.XDMFFile.__init__(self, mpicomm, filepath)
self.parameters["functions_share_mesh"] = True
self.parameters["flush_output"] = False
self.parameters["flush_output"] = True
self.path = filepath # Mimic the file path attribute from a `file` returned by `open`
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment