Error while simulating electrostatic problem for a multi-dielectric system using Gridap

I am trying to solve an electrostatic poissons equation with 2 different materials. I am trying to compute the electric potential for a sphere (ϵr= 10) embedded inside a box (ϵr= 1).

I am getting the below error at the affineFEoperator step:

julia> op = AffineFEOperator(a,b,U,V)
ERROR: promotion of types VectorValue{3, Float64} and Int64 failed to change any arguments

I am using the below code:

using Gridap, GridapGmsh, Gridap.Fields, Gridap.Geometry
const epsr_box = 1.0 ;  # Relative electric permittivity for box
const epsr_sphere = 10.0 ;  # Relative electric permittivity for sphere

model = GmshDiscreteModel("box_sph.msh");  # geometry contains a sphere domain and a box excluding the sphere i.e.. box-sphere
order = 1;
reffe = ReferenceFE(lagrangian,Float64,order);
V = TestFESpace(model,reffe,dirichlet_tags = ["HVelec", "GNDelec"]);
U = TrialFESpace(V, [2.0, 0.0]); # apply voltage at the two opposite faces of the box

degree = 1;
Ω = Triangulation(model);
dΩ = Measure(Ω,degree);

labels = get_face_labeling(model);
dimension = num_cell_dims(model);
tags = get_face_tag(labels,dimension); # get all the tags in the geometry
const sphere_tag = get_tag_from_name(labels,"sphere"); # get tag of sphere geometry
const box_tag = get_tag_from_name(labels,"box"); # get tag of box geometry

function epsr(tag)
    if tag == box_tag
        return epsr_box
    elseif tag == sphere_tag
        return epsr_sphere
    end
end

a(u,v) = ∫(  (∇(v)).(epsr*(∇(u))) )dΩ
f(x)=0.0; # assume laplace equation for now and later to have f as f(x,y,z)
b(v) = ∫(v*f)dΩ ;
op = AffineFEOperator(a,b,U,V)
uh = solve(op)

Welcome Balaji, Please read this topic to make it easier to help you. Please edit you post and write your julia code within quoted code.

(I edited the post to quote the code.)

This doesn’t look right. The dot product operator is (unicode dot operator, which can be typed with \cdot<TAB> at the julia> prompt or in editors with appropriate plugins), not . (period). A dot product can also be written dot(∇(v), epsr*(∇(u))) without the Unicode operator, but Gridap’s API is pretty Unicode-heavy in general.

See also the Gridap Poisson tutorial.

Also, I think you need to tell it to pass the tags to epsr more explicitly here. See e.g. how the tags are passed in the Gridap linear elasticity example. I think that is the source of your error message — from the way you are using it, Gridap thinks that epsr must be a function epsr(x) of the position vector x, which causes an error because epsr expects an integer tag argument.

Thank you very much. I am relatively new to julia programming. I modified the original code using the dot() operator instead of “.” as you suggested .I am writing the modified section of the new code below.

function epsr(tag)
    if tag == box_tag
        return epsr_box
    elseif tag == sphere_tag
        return epsr_sphere
    end
end
a(u,v) = ∫( dot(∇(v) , (epsr∘(∇(u),tags))) )*dΩ

I am getting a different error but at the same step as before:

julia> op = AffineFEOperator(a,b,U,V)
ERROR: MethodError: no method matching epsr(::VectorValue{3, Float64}, ::Int8)
Closest candidates are:
  epsr(::Any)

Also, I went through the documentation on multi-material elasticity modeling: 3 Linear elasticity · Gridap tutorials. I am still unclear on the syntax part. because, in the example problem the function σ_bimat has arguments (ε,tag) of which ε is itself an operator defined by 1/2(∇u+(∇u) '). Whereas in my case, the function epsr has only tag as the argument, so how can I write the weak form.
I request you to help me with 2 things.

  1. The correct modification of the code.
  2. Any link or documentation or youtube material that will help me learn the syntax of these operators/functions (as I am a beginner).

That’s because your epsr function should now take ∇(u) as the first argument since you are composing it with (∇(u),tags), e.g.

function epsr(∇u, tag)
    if tag == box_tag
        return epsr_box * ∇u
    elseif tag == sphere_tag
        return epsr_sphere * ∇u
    end
end

(similar to the tutorial I linked),