gmsh mesh and CellState projections

[christopherfear]

I am looking to project a CellState from the triangulation of a quad mesh created in gmsh.

Triangles created in gmsh/ a domain of quads made in Gridap (with the domain and partition functions) work with no issues.

The MWE of the code is as follows:

`using Gridap
using LinearAlgebra
using GridapGmsh

model = GmshDiscreteModel("C:/Users/ttcf5/gui.msh")

order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe)
U = TrialFESpace(V)

degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

f = CellState(1.0, dΩ)

function project(q,model,dΩ,order)
reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe,conformity=:L2)
a(u,v) = ∫( u*v )dΩ
l(v) = ∫( vq )*dΩ
op = AffineFEOperator(a,l,V,V)
qh = solve(op)
qh
end

fh = project(f, model, dΩ, order)`

The error is:

ERROR:
It is not possible to evaluate the given CellState on the given CellPoint.
a CellState can only be evaluated at the CellPoint it was created from.
If you want to evaluate at another location, you would need first to project the CellState
to a FESpace (e.g. via a L2 projection).

Am I misunderstanding something?
Could anyone help me to solve this problem?

[kishore-nori]

Hi,

What you have should work (just to make sure, I ran and checked, but on a different model, but that shouldn't matter). However, it might have so happened you re-ran dΩ = Measure(Ω,degree) but did not re-run the f = CellState(1.0, dΩ) once dΩ was updated and this causes a problem, as although the new dΩ is effectively the same but not identical (i.e. it is not the same object anymore, and doesn't satisfy ===, hence it results in the error you bumped into). You can restart the julia environment and re-run everything sequentially, and it should work. Having said that, a safer approach would be to have f = CellState(1.0, dΩ) inside your project function (better to also have degree = 2*order and dΩ = Measure(Ω,degree)), which will ensure f to be strictly consistent with dΩ and avoids the earlier discrepancy.

[christopherfear]

Hi,

Thanks for your suggestions. I attempted to place the suggested lines into the project function but was left with the same error output even after restarting the Julia environment. Once again it worked for a triangle mesh/quads generated in Gridap, was your test model definitely a quad mesh imported from Gmsh?

[kishore-nori]

Hi,

Yes.. I had tried it with the Gmsh generated mesh (demo.msh) taken from the Gridap Tutorials repository, and it worked, but this is not a Quad mesh but a Tet. Please check if the below works for your Gmsh model. If the issue still persists, you can attach the a minimal Gmsh file where you have this error here.

using Gridap
using LinearAlgebra
using GridapGmsh

model = GmshDiscreteModel("demo.msh")

order = 1
Ω = Triangulation(model)

function project(order,Ω)
  degree = 2*order
  dΩ = Measure(Ω,degree)
  f = CellState(1.0, dΩ)
  reffe = ReferenceFE(lagrangian,Float64,order)
  V = FESpace(model,reffe,conformity=:L2)
  a(u,v) = ∫( u*v )dΩ
  l(v) = ∫( v*f )*dΩ
  op = AffineFEOperator(a,l,V,V)
  fh = solve(op)
  fh
end

fh = project(order,Ω)

[kishore-nori]

Yeah you are right, I just tried it with a 2D quad mesh generated from Gmsh and I can reproduce the error. I think.. I understand the reason for this error, I ll try to fix it.

[christopherfear]

Thank you very much for the bug fix, that has allowed to initial projection to run.

However I think the bug goes further. I am incrementally changing the state in my full solver code. After I use the update_state! function (which takes a function comparing a CellState and OperationCellField and outputs the larger of the two) the same error returns upon attempting another projection of the result as the code loops round, as I assume the 'updated' state isn't compatible with the initial measure. I assume this is part of the same bug and wondered if you could offer further assistance?

[kishore-nori]

Thank you for reporting, I can certainly look into it and update the PR so that it can handle it. Can you provide an MWE of your issue here, so that I can explore where the issue is coming from?

[christopherfear]

The MWE I created is:

using Gridap
using LinearAlgebra
using GridapGmsh

model = GmshDiscreteModel("C:/Users/ttcf5/script_ninesegments.msh")

order = 1
Ω = Triangulation(model)
count = 1

function new_EnergyState(f,y)
  if y >= f
    f = y
  else
    f = f
  end
  true,f
end

function project(order,Ω)
  degree = 2*order
  dΩ = Measure(Ω,degree)
  if count == 1
    global f = CellState(0.0,dΩ)
  else
    global f = f
  end
  reffe = ReferenceFE(lagrangian,Float64,order)
  V = FESpace(model,reffe,conformity=:L2)
  a(u,v) = ∫( u*v )dΩ
  l(v) = ∫( v*f )*dΩ
  op = AffineFEOperator(a,l,V,V)
  fh = solve(op)
  fh
end

fh = project(order,Ω)
y = CellField(0.5, Ω) ⊙ 1 #make an example OperationCellField

update_state!(new_EnergyState,f,y)
count = count .+ 1
fh = project(order,Ω)

My full code iterates through multiple times but for the MWE I just manually typed the update_state! and then a second projection, the error of

It is not possible to evaluate the given CellState on the given CellPoint.
a CellState can only be evaluated at the CellPoint it was created from.
If you want to evaluate at another location, you would need first to project the CellState
to a FESpace (e.g. via a L2 projection).

Appears when trying to compute the second projection with the updated cellstate

[kishore-nori]

Thanks a lot for providing the MWE.

The problem here is not exactly related to the the earlier bug actually. The situation is that quadrature is not an internal state that is stored inside Gridap but generated on the fly, when user runs Measure(...), this creates the required setup for quadrature. But when we re-run Measure(...) with exactly identical arguments like done above, it doesn't result in the identical object as the vectors storing the quadrature points and weights are new vectors although holding the same values and hence the check for compatibility with your global CellState and quadrature is not met (as the compatibility checks if the CellPoints are julia identical (===), which doesn't check if the quadrature is effectively the same.., checking this is more expensive). This is natural, for example,

struct StoreVec
 x
end 

a = StoreVec([1,2])
b = StoreVec([1,2])

a === b # false 
a == b # is also false

(Measure also creates such vectors for weights and quadrature points for the reference cell).
Of course the above are essentially the same thing, but to check that, we need to check through all the values of the vectors inside and this gets more expensive if there are multiple attributes inside the struct. So, if the order and degree can decided and fixed for all your updates, we can write the function as follows, which fixes it,

order = 1
degree = 2*order
dΩ = Measure(Ω,degree)
f = CellState(0.0,dΩ)
reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe,conformity=:L2)

function project(f, V, dΩ)
  a(u,v) = ∫( u*v )*dΩ
  l(v) = ∫( v*f )*dΩ
  op = AffineFEOperator(a,l,V,V)
  fh = solve(op)
  fh
end

fh = project(f, V, dΩ)
y = CellField(0.5, Ω) ⊙ 1 #make an example OperationCellField

update_state!(new_EnergyState,f,y)
count = count .+ 1
fh = project(f,V,dΩ)

The above can be made more nicer if V and dΩ are abstracted by writing project as a struct, all being said if your order remains the same, if in case it changes, have a large enough degree to begin with to keep using the above.

[christopherfear]

Ah okay, very interesting. Worked perfectly and has allowed by MWE and full solver to run. Thanks very much for all the help with the original bug and my issues.

[kishore-nori]

That's great to hear. Thank you reporting and also for the feedback on bug issue being resolved after the PR. The PR will be pushed to master soon after the relevant tests and reviews are complete and I ll close this issue.