Hi everyone,
Im trying to solve an electromechanics problem in which I apply a potential in an internal facet that lies in the middle of my mesh, while leaving the left side fixed in the mechanical space.
So far I’ve managed to solve the electromechanics problem placing the potential at the top face, but when trying to apply the boundary condition on the inner facet, I get the following error in the variational formulation: Discontinuous type Argument must be restricted.
Here is the boundary condition set up:
import numpy as np
import ufl
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.io import gmshio
import dolfinx.io
from dolfinx import *
import dolfinx.fem.petsc
import time # To measure time taken
#import matplotlib.pyplot as plt
import dolfinx
from ufl import indices
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical, set_bc)
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_tags, facet_tags = gmshio.read_from_msh('MeshingGMSH/MyTest2.msh',mesh_comm, gdim=3)
# Function spaces
V_me = ufl.VectorElement("CG", domain.ufl_cell(), 1, 3) # Vector elems for the trial function of the mechanical problem
V_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # Scalar elems for the trial function of the electro problem
W_elem = ufl.MixedElement(V_me,V_el) # Mixed elements
W = fem.FunctionSpace(domain, W_elem) # Mixed function space
# Test and trial functions
test_fun = ufl.TestFunction(W) # Test function space
(vx, vphi) = ufl.split(test_fun)
trial_fun = ufl.TrialFunction(W)
(dx, dphi) = ufl.split(trial_fun)
x_function = fem.Function(W)
(u,phi) = ufl.split(x_function) #Llamarlos u y phi
# Mechanical subspace
W0, _ = W.sub(0).collapse() # Create a subspace from W, which corresponds to the mechanical function space (W(0))
u_bc = Function(W0) #Use a function in W0 to apply the boundary condition
u_bc.x.array[:]=0 # Set the value for the boundary condition
# Electrical subspace
W1, _ = W.sub(1).collapse() # Create a subspace from W, which corresponds to the electrical function space (W(1))
phi1_bc = Function(W1) #Use a function in W1 to apply the boundary condition
phi2_bc = Function(W1) #Use a function in W1 to apply the boundary condition
phi1_bc.x.array[:]= 0# Set the value for the boundary condition
phi2_bc.x.array[:]= 0.05# Set the value for the boundary condition
fdim = domain.topology.dim - 1
facets_left = facet_tags.find(2)
dofs_left = fem.locate_dofs_topological((W.sub(0),W0),fdim, facets_left) # Find the degrees of freedom (dofs) on the corresponding facets
#facets_top = mesh.locate_entities_boundary(domain,fdim, top) # Locate facets on the right side of the beam
facets_top = facet_tags.find(3)
dofs_top = fem.locate_dofs_topological((W.sub(1),W1),fdim, facets_top) # Find the degrees of freedom (dofs) on the corresponding facets
#facets_mid = mesh.locate_entities_boundary(domain,fdim, mid) # Locate facets on the right side of the beam
facets_mid = facet_tags.find(1)
#dofs_mid = fem.locate_dofs_geometrical((W.sub(1),W1), mid) # Find the degrees of freedom (dofs) on the corresponding facets
dofs_mid = fem.locate_dofs_topological((W.sub(1),W1),fdim, facets_mid) # Find the degrees of freedom (dofs) on the corresponding facets
facets_bottom = facet_tags.find(4)
dofs_bottom = fem.locate_dofs_topological((W.sub(1),W1),fdim, facets_bottom)
bc_mech = fem.dirichletbc(u_bc, dofs_left, W.sub(0))
bc_ele1 = fem.dirichletbc(phi1_bc, dofs_top, W.sub(1))
bc_ele2 = fem.dirichletbc(phi2_bc, dofs_mid, W.sub(1))
bcs = [bc_mech, bc_ele1, bc_ele2]
metadata = {"quadrature_degree": 4}
dV = ufl.Measure('dx', domain=domain, metadata=metadata)
ds = ufl.Measure("dS", domain=domain, subdomain_data=facet_tags, subdomain_id=1 )
Here is the mesh file with markers:
I can also post the variational formulation, but I believe Im making some mistakes while defining the boundary conditions.
Thank you very much in advance