Deleting elements in tensor mechanics

[ngrilli]

Dear Moose Users,

I was looking at the possibility of deleting elements, when a certain condition, for instance temperature increase, is satisfied.
I have already implemented models where the element stiffness is degraded, however, I would like to compare them
with a model where certain elements are actually removed from the finite element problem.

I have seen something similar is done in XFEM by creating the EFA mesh.
My understanding is that in the XFEM module, a copy of the original mesh is created and manipulated by creating cuts,
then the stress divergence kernel is calculated on the manipulated mesh. But I may be wrong here.

Can someone point me to a strategy to try the element deletion?
Is the creation of a separate mesh necessary?
I would be happy if I could do the deletion just based on element index.

Thank you very much in advance.

Best Regards,
Nicolò Grilli
National University of Singapore

[MatthewMavis]

Hi all,
I'm just jumping on this thread to say I'm also interested in this as I'm looking into topology optimisation using moose. I have developed some code to do this based on element IDs though it is currently very rough (https://github.com/MatthewMavis/moose/blob/TopOpt/framework/src/meshgenerators/ElementDeleter.C).

You should be able to hook into the ElementDeletionGeneratorBase.C like I have which will delete the elements, though one problem I'm currently having with this is that the element ID's are reassigned meaning any data from the simulation cannot be used with the original mesh. So if anyone has any ideas on how to preserve the original IDs it would be a great help.

I hope this helps in some way.

Best Regards,
Matt Mavis

[ngrilli]

Dear Matthew,
Thank you very much for your reply.
I will look into the details of your code and give you feedback.
For sure the issue with the elements ID is a difficult one.
Nicolò

[WilkAndy]

I think this would be useful to me too, but don't have time to implement it. I wrote the BlockDeleter, which may give you some hints about what to do, if you choose to do it...

[ngrilli]

Dear Andy,
Thank you very much for your reply.
Indeed, I think if block with specific IDs are deleted, the other blocks will maintain the original ID.
So this may be a better solution for me.
I will let you know how it goes
Nicolò

[ngrilli]

Dear Moose Users,

I developed the following strategy to delete elements:

1. Generate a new block with new id using ElementSubdomainIDGenerator
2. Deleting that block using BlockDeletionGenerator
3. Of course this will require a restart of the simulation, I restart the simulations at the point where I want to delete the elements using the
   SolutionUserObject and SolutionFunction as explained here:
   https://mooseframework.inl.gov/modules/porous_flow/restart.html
4. I have tested and seems the variables are assigned correctly at the correct spatial position in the elements remaining after restart
   Of course the Element Global ID changes

You can find a simple example here:
https://github.com/ngrilli/c_pfor_am/tree/main/test/tests/DeleteElementsRestart

Now the question is: will such a method preserve the Material Properties declared in the material blocks?
For instance, I have plastic deformation and I don't want it to be reinitialized to zero after restart.
Can you envisage a method to do so?
I mean restarting the simulation with elements deleted, but preserving the Material Properties at the correct spatial positions in the remaining elements.

Best Regards,
Nicolò Grilli
National University of Singapore

[jiangwen84]

@ngrilli The element deletion method is somehow on our radar to work on in the future. Instead, we have started to implement element activation for AM application recently, see #16008 It does the opposite thing to deletion, but the method would be similar. It handles material initialization as you asked about in that PR.

[jiangwen84]

@dschwen and @dewenyushu they can comment on this better.

[ngrilli]

Dear Jiang Wen,

Thank you for your answer, it would be great to have that feature in the future, but I guess it's complicated to implement.

For the moment, I am testing the restart with SolutionUserObject and SolutionFunction.
I am trying to understand if the material properties in the non-deleted elements, like plastic deformation gradient,
are initialized from the values of the previous solution or are reinitialized by the Material Class (like for a new simulation).

I am sure the plastic deformation gradient is taken from the previous solution with the "basic restart" method,
when the mesh before and after coincide,
but not sure if that is still valid with SolutionUserObject and SolutionFunction.

I will let you know what I find.

Best Regards,
Nicolò Grilli

[dewenyushu]

Our element activation capability should be readily extendable to element deletion. This is because we have "active" and "inactive" subdomains, during element activation, we move elements from "inactive" ->"active" subdomain. Element deletion will be the opposite (i.e., move elements from "active" ->"inactive" subdomain), which seems to be a minor change to #16008 .

The current PR already takes care of material initialization, updates boundary, does not require any change to the existing physics modules, and is computationally efficient since only the "active" dofs will appear in the system.

[WilkAndy]

Thanks @dewenyushu , this is a nice addition to MOOSE. In the past i have simply set the solid-mechanical stiffness of "deleted" elements to zero, or the thermal conductivity to zero, or the permeability to zero, etc (depending on the physics). But to be able to actually delete the elements is much nicer.

[ngrilli]

Dear Dewen,

Thank you very much for your answer,
I will try your modifications and give you feedback.
Yes, I think element deletion is somehow easier than insertion.

I have tried using simulation restart with SolutionUserObject and SolutionFunction,
but the Material Properties are reinitialized in the elements that are not deleted,
therefore this is not helpful for history dependent materials like plasticity.

It would be great to have your modifications released in MOOSE.

Best Regards,
Nicolò

[dewenyushu]

Sounds great! I suppose the solution provided in the PR can help with the material property issue. Let me know how it goes.

[ngrilli]

Dear @dewenyushu

I have tested your code modification for element activation/deactivation.
It works well in terms of moving elements from the "inactive block" to the "active block".

However, I see that the stress is still calculated on both active and inactive domains
and the boundary conditions I impose on the active domain are inducing stress and displacement also on the inactive domain.

This happens also in your example:
check_element_addition_by_variable.i

You mentioned in your code only the active DOF will appear,
so perhaps I am not understanding well.
This is the test I was looking at:
https://github.com/ngrilli/c_pfor_am/tree/main/test/tests/AddElements
Basically half the domain is active, half is inactive.
I apply displacement on a boundary of the active and I see the inactive domain moving as well.
Therefore the DOF of the inactive seem not be inactive.

Thank you very much in advance for your answer.

Best Regards,
Nicolò

[dewenyushu]

@ngrilli thanks for trying this!

You will need to make sure your variables (displacements and stresses), kernels (TensorMechanics), and materials are block-restricted, e.g., block = '1'.

MOOSE will throw out error about no kernel is defined in the inactive block, so you will need to add the following in your input as well.

[Problem]
  kernel_coverage_check = false
[]

[ngrilli]

Dear @dewenyushu

Yes, your suggestion worked and your code seems to add elements fine for my application.
I have created some test cases where I see the stress appearing in new elements under load:
https://github.com/ngrilli/c_pfor_am/tree/main/test/tests/AddElements

I will now try to modify the code to move elements from the active to the inactive block.
I am still studying the part of the code that updates the boundary sets and not sure that has to be changed.
For the moment I will modify only the function
ActivateElementsUserObjectBase::execute()

I will keep you updated about the results.
Thank you very much again for your help.

Best Regards,
Nicolò

[dewenyushu]

Sounds great!

Before looking deeper/trying to change the code, I would suggest trying the following:

• Define your kernels, variables, materials on the block that you want to delete elements from, say block = '2'
• Then in the ActivateElementsCoupled user object, set active_subdomain_id = '1'

It may just work. The element activation code is basically moving elements from one subdomain to another. So in theory it should not prevent us from doing element deletion without much code change. Outside of the user object, we can always define our FEProblem in the subdomain that the elements are moved to or from by using the block restricted capability.

[ngrilli]

@dewenyushu @MatthewMavis @jiangwen84

Dear Dewen,

I have completed and tested an implementation where I can do both element activation and deactivation based on a condition:
https://github.com/ngrilli/c_pfor_am

It is a minor modification of your code.
As you suggested I move elements from the inactive block to the active block and vice versa
based on the condition.

In this way the plastic deformation is reinitialised properly when elements are reactivated.

The "moving interface" is not yet updated correctly during element elimination.
I think more modifications are needed for this.
However, I am not applying any boundary condition on the "moving interface"
and my understanding is that there is stress free boundary condition at the boundary between the
active block and deactivated block.

I am happy if someone can check and test my implementation
I think it is quite interesting for the mechanics community.

Best Regards,
Nicolò Grilli
National University of Singapore

[dewenyushu]

Thanks @ngrilli for sharing your branch with us.

I looked at your example c_pfor_am/test/tests/ElementAddDelete/ElementAddDeleteDispX.i. One issue may result in incorrect moving_interface update is that you may have specified the wrong Mesh/sidesets/normal. It seems that you have left and right subdomains. Therefore, the initial moving_interface should have normal = '1 0 0' (+x direction) instead of normal='0 0 1'(+z direction). Hope this would fix the issue.

Nevertheless, I am happy that most of the part works for your application!

[Traiwit]

@ngrilli @dewenyushu

Hi guys,

I am very new to MOOSE, been using it for less than a month, but I found it very powerful.

I've been working on mine stimulation, right now I can do the blockdeletion for a simple model
However, I would like to know if it is possible to add the block/element back to the void space (backfill) without altering the strain on surrounding area (strain-free).

to make it clearer, here's the workflow:
create a solid geometry -> remove a block -> run a equalibirum simulation with gravity activate -> the void is deformed due to stress
-> backfill the material back into the void without altering the deformation from the previous step.

Another question I have is: when I use 'BlockDeletionGenerator', the block is gone, the issue would also be that: how can i locate the position of the backfill if I no longer have the deleted block in the system.

I read @ngrilli reply back in Jan, but i am not sure if it will works for my problem, could you please clarify, thanks!.

Thank you

Kind regards,
Trai

below is my script for blockdeletion

[Mesh]
[./gen]
type = FileMeshGenerator
file = equlstep_step0.e
# show_info = true
use_for_exodus_restart = true
# exodus_extra_element_integers = 's11_aux s22_aux s33_aux'

[]
[mergedBox]
type = RenameBlockGenerator
input = gen
old_block_id = '2'
new_block_id = '1'
[]
[ed0]
type = BlockDeletionGenerator
input = mergedBox
block_id = 1
[]
[]

[jlinBE]

Hi,

I'm working on the same problem alongside Traiwit, and I think all we need to do at this point is: 1. let the pre-backfill simulation reach equilibrium. 2. import the final total strain in a new simulation, as the (negative) eigenstrain.

I am pretty sure this would be possible, but I don't know how to make the Moose syntax work this way. Does anyone have an example to share?

[WilkAndy]

I suggest starting a new discussions thread. It's likely people have "tuned out" to this discussion, and it's not about deleting elements anymore.

I think what you want is:

1. At the start of each timestep, set the Young's modulus of certain elements to "zero" (as you excavate the mine). The "certain elements" are defined by a time-dependent Function, or AuxVariable, or whatever. The "zero" is a small number: it should not be exactly zero, otherwise the displacements won't be uniquely defined. Different elements will be "set to zero" at each timestep. You've already achieved this (have you looked at https://mooseframework.inl.gov/modules/porous_flow/coal_mining.html, by the way).
2. At the start of each timestep, set the strain in certain elements to zero. The "certain elements" are defined by another time-dependent Function, etc. Zeroing the strain should not impact the positions of the nodes in any element. Its purpose is to ensure the stress in these elements is zero (at the start of the timestep). This means an alternative would be to introduce a new reference strain, strain_ref, which is the strain in those elements at the start of the timestep, and calculate stress using stress = E * (strain - strain_ref). You haven't be able to achieve this, since there appears to be nothing in MOOSE that enables either of these ideas.

[jlinBE]

Hi Andy,

Thanks for taking the time to reply to us! Trai is writing up a new thread.

This is exactly right! And in fact we did it while referencing the mech-only version of coal_mining in tensor_mechanics/examples/coal_mining :-)

The key differences in our case is we are solving for steady-state at every step instead of transient, and we're using the base tensormechanics stress/strain classes instead of cosserat.

Since we are solving for steady-state, we can't define time-dependent functions. Other than that, this is exactly what we're trying to do. Instead of taking the strain at the start of the timestep, we're trying to take it from the steady-state solution of the previous simulation.

I think the eigenstrain-related materials (ComputeEigenstrain and ComputeVariableEigenstrain) can probably do what I want, but not sure about the syntax.

[WilkAndy]

OK, sounds good. Instead of doing steady-state, you could do a transient quasi-static simulation, which avoids having to read/write files.

[ngrilli]

Dear @Traiwit

You mention that you want to reduce to zero the stiffness of some of your elements
but it seems in the input file blocks attached to your messages that you do not have any object to do so.

You will need to make a new object that degrades the stiffness based on some criterion,
for instance the value of a coupled variable.

You can see an example code for stiffness degradation in my repository:
https://github.com/ngrilli/c_pfor_am/blob/main/src/materials/ComputeElasticityTensorMelting.C
In my case the coupled variable is the temperature.

Best Regards,
Nicolò Grilli
National University of Singapore