Cohesive zone modeling, need help

[Andrea Rovinelli]

Dear all,
I need to setup a cohesive zone model to study intra-granular cracking. In this model cracking can occur only at grain boundaries and to implement this behavior I was thinking to use cohesive elements. I would have XFEM because I don't need an arbitrary crack path, but just grain boundary cracking.
First I'm not sure if MOOSE can handle cohesive elements, i tried to look into libmesh supported element types (https://libmesh.github.io/doxygen/classlibMesh_1_1Elem.html) and I don't see specific elements allowed to have zero thickness. In fact when I try a simple 2D case in which I add some zero thickness QUAD4 elements moose complains saying that the element has zero volume.

I know someone from the MOOSE development team is working on this feature and that at some point it will be released, but we need to work on this project.

Any hint and help would be really appreciated.

Thanks in advance


Andrea (ANL)

[Wen Jiang]

Hi Andrea,

I spent a few days last months to develop CZM in MOOSE. I found it was NOT trivial work and the implementation is far from completion. Anyway here is my branch if you would like to take a look at what I have done (https://github.com/jiangwen84/moose/tree/czm_moose)

Basically, what I did was

1) write a new mesh class: it reads from a standard mesh file and at cohesive interface it adds nodes and changes element connectivities.

2) use interface kernel to implement the cohesive zone law. 

3) use boundary material system to store stateful properties. 

I am not sure if I have time to work on this soon, but I can help you if you would like to continue my work.

Regards,

Wen

[Andrea Rovinelli]

Dear Wen,
thanks for your reply. Our idea is to start from you work and develop this capability.
First let me say that I checkout your czm_moose branch, and after recompiling libmesh and running test I had the follwoing errors

materials/stateful_prop.spatial_bnd_only.................................................... FAILED (EXODIFF)
materials/stateful_prop.stateful_copy....................................................... FAILED (EXODIFF)
materials/stateful_internal_side_uo.test.................................................... FAILED (EXODIFF)
materials/stateful_prop.spatial_test........................................................ FAILED (EXODIFF)
userobjects/internal_side_user_object.get_neighbor_test..................................... FAILED (EXODIFF)

I don't know if this relevant or not to what you were doing.

Also I noticed you have a test input file czm.i  but the related czm0.e file is missing. If you have it can you upload it?

Then I have a more fundamental question, why you are inserting additional nodes in between the interfaces rather then implementing cohesive elements in libmesh, it shouldn't be this the most natural way to do it? There is a practical reason for this choice?

Thanks again

Andrea

[Wen Jiang]

On Friday, February 2, 2018 at 5:32:45 PM UTC-7, Andrea Rovinelli wrote:

Dear Wen,
thanks for your reply. Our idea is to start from you work and develop this capability.
First let me say that I checkout your czm_moose branch, and after recompiling libmesh and running test I had the follwoing errors

materials/stateful_prop.spatial_bnd_only.................................................... FAILED (EXODIFF)
materials/stateful_prop.stateful_copy....................................................... FAILED (EXODIFF)
materials/stateful_internal_side_uo.test.................................................... FAILED (EXODIFF)
materials/stateful_prop.spatial_test........................................................ FAILED (EXODIFF)
userobjects/internal_side_user_object.get_neighbor_test..................................... FAILED (EXODIFF)

I don't know if this relevant or not to what you were doing.

Like what I said, my code was not complete yet and it broke some existing tests. Basically, MOOSE would calculate boundary materials only if there exists integrals bc. To use the boundary materials on the cohesive interface, I made some hacky codes in the framework which need to be corrected. 

 

Also I noticed you have a test input file czm.i  but the related czm0.e file is missing. If you have it can you upload it?

I will send it to you shortly. 

 

Then I have a more fundamental question, why you are inserting additional nodes in between the interfaces rather then implementing cohesive elements in libmesh, it shouldn't be this the most natural way to do it? There is a practical reason for this choice?

There are two ways that you can use to implement CZM. One is to insert cohesive element with zero thickness. The other one is to treat the cohesive zone as an interface. The second one is more like a Discontinuous Galerkin type way. You add additional nodes and then constraint them. MOOSE implements the interfaceKernel based on DG and it should be a good fit for CZM. Regarding the first approach, we might need to add special element type in MOOSE and modify its integration/assembly staff. It seems to have more work to do in my opinion. 

Let us talk offline for the implementation details.

Regards,

Wen     

[Andrea Rovinelli]

Dear Wen,
I worked on the base you provided me and I have been able to implement the Xu Needleman model for the cohesive zone.
Code is here https://github.com/arovinelli/moose/tree/czm_wen .
You can compile it and execute the test moose/modules/tensor_mechanics/test/tests/cohesive_zone/czmXN.i

Now that a non-trivial traction separation law works we should find a better way to implement this interface (ie. starting from an already disjoint mesh with hacking during runtime) to make it work with large deformation and convoluted interfaces, like in the case of grain boundaries.
In fact, the current implementation does not work with triple points (see the test czm_junction_2D.i ).
Unfortunately this i beyond my Libmesh and Mosse knowledge. I know Derek was working on it, but I a haven't heard from him in a while.


Andrea

[Wen Jiang]

Hi Andrea,

I believe Derek has been working on this and it should be fairly easy to implement CZM for grain boundary based on his implementation.

Since we have decided to use FaceFaceConstraint for czm, I do not want you to waste your time to develop czm based on my old prototype implementation using interfaceKernel :-)

Derek: do you have any progress to share with us? Thanks.

Regards,

Wen

[Gaston, Derek R]

I'm here - but many things have gotten in the way (as usual).  Let me review this and get back to you soon.

Derek

--
You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/8945f6bb-982e-4a5d-aa10-611dfdf1b289%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Andrea Rovinelli]

Wen thanks for the reply.
Derek thanks also to you, I'll wait for your update.

All the Best
Andrea

On Wednesday, April 4, 2018 at 2:11:55 PM UTC-5, Derek Gaston wrote:

I'm here - but many things have gotten in the way (as usual).  Let me review this and get back to you soon.

Derek

On Wed, Apr 4, 2018 at 11:39 AM, Andrea Rovinelli <andrea.r...@gmail.com> wrote:

Dear Wen,
I worked on the base you provided me and I have been able to implement the Xu Needleman model for the cohesive zone.
Code is here https://github.com/arovinelli/moose/tree/czm_wen .
You can compile it and execute the test moose/modules/tensor_mechanics/test/tests/cohesive_zone/czmXN.i

Now that a non-trivial traction separation law works we should find a better way to implement this interface (ie. starting from an already disjoint mesh with hacking during runtime) to make it work with large deformation and convoluted interfaces, like in the case of grain boundaries.
In fact, the current implementation does not work with triple points (see the test czm_junction_2D.i ).
Unfortunately this i beyond my Libmesh and Mosse knowledge. I know Derek was working on it, but I a haven't heard from him in a while.


Andrea

--
You received this message because you are subscribed to the Google Groups "moose-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

[Andrea Rovinelli]

Dear Derek,
did you had a chance to look at the constraint?

Thanks
Andrea

[Andrea Rovinelli]

Wen i've been able to implement a mix mode cohesive law.
Now I need to couple additional displacement components to properly impose the offidagonal jacobian.
Could you tell me which commands I need to insert in the .h and .C file.

Also, Derek have you been able to recheck the new constraint?

All the best

On Wednesday, April 4, 2018 at 2:33:46 PM UTC-5, Andrea Rovinelli wrote:

[Wen Jiang]

Are you using InterfaceKernel? If so, you can use 

InterfaceKernel::computeElementOffDiagJacobian(unsigned int jvar)

and 

InterfaceKernel::computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int jvar)

to provide off-diagonal jacobian. 

[Alexander Lindsay]

Unless you're doing something non-standard, the interfaces should be computeQpJacobian(Moose::DGJacobianType type) and computeQpOffDiagJacobian(Moose::DGJacobianType type, unsigned int jvar)

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/53084002-6a88-433b-bfac-098388144cf7%40googlegroups.com.

[Andrea Rovinelli]

Wen and Alex, thanks for the help.
Yes I'm using the interface kernel. I have one question, if on both side of the interface I couple the same variable (let's say displacement_x), will I receive two different variable's identifier when calling

  _disp_1_var(coupled("disp_1")),
  _disp_1_neighbor_var(coupled("disp_1_neighbor")),

or they would be the same?

thanks

[Andrea Rovinelli]

I also have another question, which is more related to efficency.
When computing residual, jacobian and offdiagjacobian, do I have to recompute all the intermediate variables required to obtain those values or mosse it smart enogh.
Specifically I have a few method used to compute the residual and the complete jacobian matrix. do I have to call them separetly in computeQpResidual, computeQpJacobian and computeQpOffDiagJacobian or i can compute them ones and store the result for the current iteration somewhere else inside the kernel?

Thanks

[Cody Permann]

This might be a case of premature optimization...

I would suggest that you don't worry about doing extra work at this stage. Let's get something up and running first and see if there are speed issues. You should note that of these methods that you mentioning here are scalable! As you run on more processors each process has fewer elements to work on so even if you are doing a little bit of rework, it's scalable. Any kind of caching that you implement here is going to be trading off time for space, and we don't even know what percentage of the calculation is going into the overlapping part of these routines without profiling first so it's just premature to worry about this right now.

--
You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/56beed33-c873-47af-9933-0af2c7fb0650%40googlegroups.com.

[Andrea Rovinelli]

Cody,
thanks for the suggestion. For the moment, just to be sure I'll call all the required method in each computeQP routine.
All the best.

[Alexander Lindsay]

Are you declaring two different variables in your input file? Or are the names you're passing into `variable` and `neighbor_var` the same?

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/d0fcf407-105e-4647-847f-3c7f274ee9ae%40googlegroups.com.

[Andrea Rovinelli]

Alex,
I'm declaring a single variable for each displacement component on all the grid.
the input file for the interface kernel taking care of separation looks like

  [./interface_y]
    type = InterfaceCohesiveZoneUObased
    variable = disp_y
    neighbor_var = disp_y
    disp_1 = disp_x
    disp_1_neighbor = disp_x
    disp_index = 1
    boundary = 100
  [../]

[Alexander Lindsay]

Then the IDs should be identical

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/de436d68-4746-4e63-a8bb-84b473e8707f%40googlegroups.com.

[Andrea Rovinelli]

Dear All,

I've worked on this implementation for a while and now I have a nice package.

It includes 4 components:

1. A material object (DicplacementJumpBasedCohesiveInterfaceMaterial) which reads the displacement jump across an interface (the material also convert from global to local coordinates and vice versa)
2. A user object superclass (TractionSeparationUOBase) and a derived class (CohesiveLaw_3DC) in which the Residual and Jacobian are computed in the interface coordinates
3. An interface Kernel (DisplacementJumpBasedCohesiveInterfaceKernel). The interface kernel works on 1 coordinates at a time, so we need to 3 kernels for 3D material computation
4. A mesh modifiers which split the mesh (each block becomes a region) and adds the boundary interfaces. This is based on what Wen did, I extended it to 3D and also to deal with points junctions

Everything works fine wen the interface is flat and we have only 2 blocks.

Here I add 2 examples solution, the first for pure normal mode and the second for pure shear mode (input files are czmUONormal_3D_x45_y45.i and czmUOShear_3D_x45_y45.i, respectively). Please note that this simulations are performed on rotated mesh to check that rotations and boundary normals are managed properly.

However, when the geometry of the interface becomes complex and points junctions are present the solver cannot converge.

I checked and rechecked the modified mesh and everything seems fine to me. Below there are 2 examples of exploded modified meshes with their relatives boundary interfaces (czmUONormal_3Blocks_Mesh.i and czmUONormal_3Blocks_Mesh2.i).

Shades of blue represents blocks, and shades of red different interfaces (1 for each pair of blocks). Note that having one interface for each pair of blocks is inconvenient, but I thought that having them separate may have solved the issue.

I should mention that before writing this post I tried different options to check what was wrong.

I tried with one global interface, I tried to add different materials kernels and userobject for each interface, ... The result is that everything works fine only if there are no npoint junctions.

Now I'm stuck and I would really appreciate if someone of the developers (Cody, Daniel,  Wen, Alex, Derek :) ) could review this peace of work and give me some insight on how to proceed further, or point me to the problem.

The code can be found here https://github.com/arovinelli/moose/tree/czm_wen .

Note that all this work is based on the initial model proposed by wen (all the files ending in Wen are the original ones). 

Also it would be really nice to avoid splitting the mesh inside moose but to do this I will need the constraint Derek was working. I could also try to sue some GeomSearch routine but I'm not familiar with them at all.

As a note, I know using the DG interface kernel approach might not be the best approach, but we need something working to simulate grain boundary cracking in polycrystalline aggregates. At this point in time, as we don't have alternatives we are fine with DG kernel.

Kind Regards

Andrea

[walkand...@gmail.com]

Thank you Andrea!

在 2018年5月11日星期五 UTC+2下午5:44:37，Andrea Rovinelli写道：

[Andrea Rovinelli]

Dear all,

during the weekend I spent some time in trying to figure out what's wrong with my implementation.

To make sure that the cause of the problem is not the mesh splitting algorithm, I did a step back.

I used a 2D mesh and hard-coded the splitting (i.e. I manually duplicated nodes and assigned sides to the boundary).

Mesh modifications are performed in CohesiveZoneMesh2D_JunctionTest. Tests performed with this model can be found in 

.../tensor_mechanics/test/tests/CohesiveZoneJunction.

Even in these cases the solver does not converge, therefore the problem is somewhere else.

I think the way the boundary interface is setup does not allow to properly add residuals at npoint junctions. 

Any comment or help would be really appreciated

All the best

[Wen Jiang]

I think the interfaceKernels should add residuals correctly regardless of the interface topology. I need to pull down your branch and see if I can find what causes it fail to converge. 

[walkand...@gmail.com]

Hi Wen,

is it possible that in a Material kernel, we can get some neighbor information via the coupledNeighborValue function?

I read your code and Andrea's code, for example the MaximumNormalSeparationWen, the _max_normal_separation is

calculated by disp_y and disp_y_neighbor.

I try to do this in my material kernel, but it gives me some compile error.

So how can we get the neighbor information in a material kernel correctly?

Best regards

在 2018年5月14日星期一 UTC+2下午8:36:13，Wen Jiang写道：

[Andrea Rovinelli]

Wen,
Thanks for taking a look at it,
As suggestion to make things easy, i would start from the czmTestJunction_2D.i, which uses the CohesiveZoneMesh2D_JunctionTest Mesh modifier. In this case the mesh is composed by 4 elements and is manually splitted.

let me know if you need anything else.

All the best
Andrea

[Andrea Rovinelli]

Dear Yang,
In this branch, the file framework/include/materials/Material.h
has been modified to allow the coupling of material at the interface.
You can see all the modifications we did to the original moose kernel here https://github.com/idaholab/moose/compare/devel...arovinelli:czm_wen

Probably we should create a new material class that implement this behavior, I'll do this when everything else is working, prior to merge.

Andrea

[walkand...@gmail.com]

Thank you Andrea, I will check it.

在 2018年5月14日星期一 UTC+2下午9:31:08，Andrea Rovinelli写道：

[Andrea Rovinelli]

Ok, Digging into the problem I figured out there is a sign error in the residual.

What is the appropriate sign convention for the interface kernel for Elem and Neighbor (and also for the Jacobian)
In this case my traction separation equation is written as u_neighbor-u and a positive jump (let's say in the opening direction) produces a positive traction in that direction whcih is also aligned with the normal of the interface.

On Monday, May 14, 2018 at 1:36:13 PM UTC-5, Wen Jiang wrote:

[Alexander Lindsay]

On the master side, the sign convention is the same as it would be for an IntegratedBC that falls out of your integration by parts (generally speaking). On the slave side it will be the opposite sign of what would fall out of integration by parts because the normal vector is defined with respect to the master side (e.g. pointing outward from the master volume)

--

You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/80cf94a5-ac39-4dde-b611-2d459f98d334%40googlegroups.com.

[Andrea Rovinelli]

Alex thanks for your reply.
I have one last quick question.
I would like to check my jacobian implementation utilizing the built in jacobian debugger.
I use SMP full for my preconditioner block. Now I know that BCs are discarded so I setted up ICs, but it does seems to ignore my interface kernel.
So for an interface kernel coupling the displacement jump between two blocks, which is the correct initial condition to set.
(note that the Jump, residual and Jacobian are computed in the material and then imported in the interface kernel).

Thanks

On Monday, May 14, 2018 at 5:17:42 PM UTC-5, Alexander Lindsay wrote:

On the master side, the sign convention is the same as it would be for an IntegratedBC that falls out of your integration by parts (generally speaking). On the slave side it will be the opposite sign of what would fall out of integration by parts because the normal vector is defined with respect to the master side (e.g. pointing outward from the master volume)

On Mon, May 14, 2018 at 4:00 PM, Andrea Rovinelli <andrea.r...@gmail.com> wrote:

Ok, Digging into the problem I figured out there is a sign error in the residual.

What is the appropriate sign convention for the interface kernel for Elem and Neighbor (and also for the Jacobian)
In this case my traction separation equation is written as u_neighbor-u and a positive jump (let's say in the opening direction) produces a positive traction in that direction whcih is also aligned with the normal of the interface.


On Monday, May 14, 2018 at 1:36:13 PM UTC-5, Wen Jiang wrote:

I think the interfaceKernels should add residuals correctly regardless of the interface topology. I need to pull down your branch and see if I can find what causes it fail to converge. 

On Monday, May 14, 2018 at 8:09:20 AM UTC-6, Andrea Rovinelli wrote:

Dear all,

during the weekend I spent some time in trying to figure out what's wrong with my implementation.

To make sure that the cause of the problem is not the mesh splitting algorithm, I did a step back.

I used a 2D mesh and hard-coded the splitting (i.e. I manually duplicated nodes and assigned sides to the boundary).

Mesh modifications are performed in CohesiveZoneMesh2D_JunctionTest. Tests performed with this model can be found in 

.../tensor_mechanics/test/tests/CohesiveZoneJunction.

Even in these cases the solver does not converge, therefore the problem is somewhere else.

I think the way the boundary interface is setup does not allow to properly add residuals at npoint junctions. 

Any comment or help would be really appreciated

All the best

--
You received this message because you are subscribed to the Google Groups "moose-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

[Alexander Lindsay]

What do you mean it ignores it? When I'm testing Jacobians, I always use RandomICs to make sure gradients are non-zero

[Andrea Rovinelli]

What do you mean it ignores it?

The  debugger always returns "no problem detected" even if I purposely implement a wrong Jacobian
Could this be related to the fact that the residual and jacobian are computed and stored as material property and only pulled from the kernel?
 

When I'm testing Jacobians, I always use RandomICs to make sure gradients are non-zero

I do set RanomICs on the displacement variables on the master side of the interface boundary.

[Alexander Lindsay]

On Wed, May 16, 2018 at 8:15 AM, Andrea Rovinelli <andrea.r...@gmail.com> wrote:

What do you mean it ignores it?

The  debugger always returns "no problem detected" even if I purposely implement a wrong Jacobian
Could this be related to the fact that the residual and jacobian are computed and stored as material property and only pulled from the kernel?

This makes it sound like I need to go look at your code. I don't understand what the above sentence means. Are you not still returning residuals and Jacobians from your computeQpResidual/Jacobian methods?

 

When I'm testing Jacobians, I always use RandomICs to make sure gradients are non-zero

I do set RanomICs on the displacement variables on the master side of the interface boundary.

--
You received this message because you are subscribed to the Google Groups "moose-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/4c4afbbe-4cb2-40b4-bc5d-c672b986a6f7%40googlegroups.com.

[Andrea Rovinelli]

I'm returning them from the standard methods (computeQpResidual/Jacobian ) but the coefficients needed to construct the residual and the jacobian are stored into a material property


into the cosntructor i have
_ResidualMP = &getMaterialProperty<RealVectorValue>(_residual);

and the computeQPResidual is

Real
DisplacementJumpBasedCohesiveInterfaceKernel::computeQpResidual(Moose::DGResidualType type)
{

  // weak form of the problem [test_neighbor-test]*r
  Real r = (*_ResidualMP)[_qp](_disp_index);

  switch (type)
  {
  
    case Moose::Element:
      r *= -_test[_i][_qp];
      break;

    case Moose::Neighbor:
      r *=   _test_neighbor[_i][_qp];
      break;
  }

  return r;

}





On Wednesday, May 16, 2018 at 12:46:39 PM UTC-5, Alexander Lindsay wrote:

On Wed, May 16, 2018 at 8:15 AM, Andrea Rovinelli <andrea.r...@gmail.com> wrote:

What do you mean it ignores it?

The  debugger always returns "no problem detected" even if I purposely implement a wrong Jacobian
Could this be related to the fact that the residual and jacobian are computed and stored as material property and only pulled from the kernel?

This makes it sound like I need to go look at your code. I don't understand what the above sentence means. Are you not still returning residuals and Jacobians from your computeQpResidual/Jacobian methods?

 

When I'm testing Jacobians, I always use RandomICs to make sure gradients are non-zero

I do set RanomICs on the displacement variables on the master side of the interface boundary.

--
You received this message because you are subscribed to the Google Groups "moose-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

[Andrea Rovinelli]

dear all.
While debugging the Jacobian I found a strange problem.
The material properties pulled from the interface kernel belong to the wrong element (I.e kernel is working on element 0 and material property come from element 2)

What could be the cause of the problem??
My impression is that my material object is saving properties on the wrong element, but i’m Not sure how to check.

[Wen Jiang]

In my czm branch, I used boundary material property to calculate the opening and traction force and use them in the interfaceKernel. The 'boundary' is defined on only one side of the interface. 

The material properties pulled from the interface kernel belong to the wrong element (I.e kernel is working on element 0 and material property come from element 2)

What could be the cause of the problem??

 Could you check if element 0 and element 2 are neighbor? If so, it is possible that the InterfaceKernel works on elem0 and elem2 face while fetching the material property from elem2 face.

[Andrea Rovinelli]

Wen
Thanks for hint

In my czm branch, I used boundary material property to calculate the opening and traction force and use them in the interfaceKernel. The 'boundary' is defined on only one side of the interface. 

my implementation is based on your original branch. There are more stuff but the underlying core is the same. 

 Could you check if element 0 and element 2 are neighbor? If so, it is possible that the InterfaceKernel works on elem0 and elem2 face while fetching the material property from elem2 face.

Element 2 and 0 are neighbor, however the properties coming from the MP are NOT from the neighboring face (see the mesh and kernel logs).
I can tell this by also looking at the normals of the face which are not only flipped.
I should also mentioned that I started facing convergence problems when I started using non planar interfaces where all the normals are the same.
So tehre might be something else going on here,
By the way in the image  blocks correspond to element IDs.



Mesh log
elem id = 0
>>>>>>>>>>>>>>>>>>>>>>>>>>>>

ELEM
  Elem Information
   id()=0, unique_id()=5, processor_id()=0
   type()=TRI3
   dim()=2
   n_nodes()=3
    0  Node id()=0, processor_id()=0, Point=(x,y,z)=(0.00424005, -6.30809e-18,        0)
    DoFs=
    1  Node id()=1, processor_id()=0, Point=(x,y,z)=(    -0.5,     -0.5,        0)
    DoFs=
    2  Node id()=2, processor_id()=0, Point=(x,y,z)=(     0.5,     -0.5,        0)
    DoFs=
   n_sides()=3
    neighbor(0)=2
    neighbor(1)=NULL
    neighbor(2)=1
   hmin()=0.704115, hmax()=1
   volume()=0.25
   active()=1, ancestor()=0, subactive()=0, has_children()=0
   parent()=NULL
   level()=0, p_level()=0
   refinement_flag()=DO_NOTHING
   p_refinement_flag()=DO_NOTHING
   DoFs=

>>>>>>>>>>>>>>>>>>>>>>>>>>>>NEIGHBOR_ELEM
elem id = 2
>>>>>>>>>>>>>>>>>>>>>>>>>>>>ELEM
  Elem Information
   id()=2, unique_id()=7, processor_id()=0
   type()=TRI3
   dim()=2
   n_nodes()=3
    0  Node id()=6, processor_id()=0, Point=(x,y,z)=(0.00424005, -6.30809e-18,        0)
    DoFs=
    1  Node id()=4, processor_id()=0, Point=(x,y,z)=(    -0.5,      0.5,        0)
    DoFs=
    2  Node id()=9, processor_id()=0, Point=(x,y,z)=(    -0.5,     -0.5,        0)
    DoFs=
   n_sides()=3
    neighbor(0)=3
    neighbor(1)=NULL
    neighbor(2)=0
   hmin()=0.710111, hmax()=1
   volume()=0.25212
   active()=1, ancestor()=0, subactive()=0, has_children()=0
   parent()=NULL
   level()=0, p_level()=0
   refinement_flag()=DO_NOTHING
   p_refinement_flag()=DO_NOTHING
   DoFs=

Kernel log
==============WRONG ELEMENT INFO==================================
 ELEMENT
 MP element =====> 2
 Kernel element =====> 0

==============================
 SIDE
 MP side =====> 0
 Kernel side =====> 0

==============================
 QP
 MP qp =====> 0
 Kernel qp =====> 0

==============================
 NORMALS OUT OF SYNC
full MP LocalJump
0.704115   0.710086   -0  
full Kernel LocalJump
-0.704115   0.710086   0  
============================================================
 JUMP OUT OF SYNC
full MP LocalJump
0.661071   0.0771329   0  
full Kernel LocalJump
-0.496427   0.796534   0  
============================================================
 LOCAL JUMP OUT OF SYNC
full MP LocalJump
0.520241   0   0.415107  
full Kernel LocalJump
0.915149   0   0.208346  
============================================================
 LOCAL TRACTION OUT OF SYNC
full MP Traction
42.1912   0   28.8765  
full Kernel Traction
83.7404   0   16.3529  
============================================================
 LOCAL TRACTION DERIVATIVES OUT OF SYNC
full MP TractionDerivatives
38.9081   0   -140.111  
0   69.564   0  
-28.8765   0   -26.3307  

full Kernel TractionDerivatives
7.76426   0   -139.576  
0   78.4894   0  
-16.3529   0   51.2329  

============================================================
  TRACTION OUT OF SYNC
full MP Traction
50.2122   9.62698   0  
full Kernel Traction
-47.3509   70.9772   0  
============================================================
 TRACTION DERIVATIVES OUT OF SYNC
full MP TractionDerivatives
-78.4774   87.5219   0  
-23.7124   91.0548   0  
0   0   69.564  

full Kernel TractionDerivatives
107.644   82.6865   0  
-40.5361   -48.6464   0  
0   0   78.4894  

 

[Andrea Rovinelli]

I'm wondering if it might be related to the changes

framework/src/base/ComputeMaterialsObjectThread.C

between lines 110 and 120, specially the one

_assembly[_tid]->reinit(elem, side);

[Andrea Rovinelli]

I spent a bit more time on this.
and reported a log of all the wrong element pairs between the Material and Kernel. note that even sides sometimes are mixed (qp are always ok but I can't say that this is only a lucky chance).
Now the question is: if element and side renumbering is performed, where would it be located in the code?
Note that the pulled material property are consistent with the original mesh numbering (I can tell this buy looking at element id, side and normal), so I guess that there is some pointer somewhere going crazy a pointer somewhere going crazy.

Thanks

[Andrea Rovinelli]

By the way, if someone wants to take a look at a code and pull it this is the link

https://github.com/arovinelli/moose/tree/czm_wen

Andrea

[Alexander Lindsay]

I'm losing track of what's going on here. Without boring down into the code, it seems like this customization of the InterfaceKernel system has gone a bit off the tracks. If I was to pull your branch, would I find nothing but easy-to-understand, beautiful code? I'm loathe to go down a rabbit hole for a concept that seems to be a more natural fit for an entirely different system (FaceFaceConstraints)...

--

You received this message because you are subscribed to the Google Groups "moose-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/d6e05070-1c37-4343-9e99-b9e6649bdf77%40googlegroups.com.

[Andrea Rovinelli]

Dear Alexander

I'm losing track of what's going on here. Without boring down into the code, it seems like this customization of the InterfaceKernel system has gone a bit off the tracks.

let me explain what's happening. I had issues with the jacobian and I wasn't able to understand what was wrong with math. After rechecking the equations dozens of times and making a simple FD check I come to the conclusion that the issue was somewhere else. The original model I posted on this mailing list a couple of weeks ago was divided into 3 components a material (doing all the computation on the interface) a user object containing the CZ law and an interface kernel pulling residual and jacobian coefficients from the material. To check if my jacobian was really wrong or not I did all the computation inside the interfacekernel. As soon as I did this I achieved quadratic convergence and no issues. Therefore the jacobian is right and the interface kernel behave as expected. However I need material property on the interface to keep track of state variables. Therefore I compared one by one the result of MP with the ones obtained in the interfacekernel. It comes out that when retrieving material properties something goes wrong, specifically the property pulled into the interfacekernel sometimes do not belong to the current element or side.
I don't know where this issue is originating. It might be a consequence of splitting the mesh at runtime (we discussed this in the telecon a few months ago), or it might be something else.
As far as my knowledge goes this might be related to my implementation or to the original changes that Wen did to allow a material to be neighbor coupleable, or to other changes related to the ComputeMaterialObjectThread or ComputeResidualThread. I'm not an expert of these core libraries so I can't comment on this. Those changes are what I asked someone with more experience than me to look at not all my code.
 
 

  If I was to pull your branch, would I find nothing but easy-to-understand, beautiful code?

We can discuss about beautiful code styling as much as you want. I don't think you will find my code "beautiful" as it is in the middle of my development and I'm making changes all the time to figure out what is really going on, but I don't consider myself a poor coder. Again it was not my intent to let you review all the implementation but just the changes to the basic routines which are only a few lines in 3 or 4 files (again note that I didn't do these changes).

I'm loathe to go down a rabbit hole for a concept that seems to be a more natural fit for an entirely different system (FaceFaceConstraints)...

I asked for help in converting the FaceFaceConstraint (which now is deprecated and is renamed as Mortar) from mortar to DG a number of times. I also asked for suggestion on how to do it myself as you might have more compelling things to work on (email and this discussions here) but you know better than what the answer has been until today. If you can provide me with guidance in doing the conversion from mortart to DG or penalty I can do that (this is not a moose user change, is more of adding a core functions that requires deep knowledge of the framework). HoweverI have to say that conventional Cohesive elemnents does not differe so much in their limitation to the interface kernel where the neighboorhood is considered to be the one present in the undeformed configuration (i.e. calculation are performed always  between the same faces). These is why we decided to follow the route Wen implemented months ago. Also It was our only option at that time.

I hope this clarify your question and concerns.
Kind regards

Andrea
 

[Alexander Lindsay]

On Fri, May 18, 2018 at 10:47 AM, Andrea Rovinelli <andrea.r...@gmail.com> wrote:

Dear Alexander

I'm losing track of what's going on here. Without boring down into the code, it seems like this customization of the InterfaceKernel system has gone a bit off the tracks.

let me explain what's happening. I had issues with the jacobian and I wasn't able to understand what was wrong with math. After rechecking the equations dozens of times and making a simple FD check I come to the conclusion that the issue was somewhere else. The original model I posted on this mailing list a couple of weeks ago was divided into 3 components a material (doing all the computation on the interface) a user object containing the CZ law and an interface kernel pulling residual and jacobian coefficients from the material. To check if my jacobian was really wrong or not I did all the computation inside the interfacekernel. As soon as I did this I achieved quadratic convergence and no issues. Therefore the jacobian is right and the interface kernel behave as expected. However I need material property on the interface to keep track of state variables. Therefore I compared one by one the result of MP with the ones obtained in the interfacekernel. It comes out that when retrieving material properties something goes wrong, specifically the property pulled into the interfacekernel sometimes do not belong to the current element or side.
I don't know where this issue is originating. It might be a consequence of splitting the mesh at runtime (we discussed this in the telecon a few months ago), or it might be something else.
As far as my knowledge goes this might be related to my implementation or to the original changes that Wen did to allow a material to be neighbor coupleable, or to other changes related to the ComputeMaterialObjectThread or ComputeResidualThread. I'm not an expert of these core libraries so I can't comment on this. Those changes are what I asked someone with more experience than me to look at not all my code.

This is a nice summary. Thank you. It's good to know that the InterfaceKernel, at least by itself, is doing what it should do. At this point, I would say this is a development issue. I've created a ticket on github where I suggest we move further discussion. I'm going to take a look at your branch.

Alex

 

 
 

  If I was to pull your branch, would I find nothing but easy-to-understand, beautiful code?

We can discuss about beautiful code styling as much as you want. I don't think you will find my code "beautiful" as it is in the middle of my development and I'm making changes all the time to figure out what is really going on, but I don't consider myself a poor coder. Again it was not my intent to let you review all the implementation but just the changes to the basic routines which are only a few lines in 3 or 4 files (again note that I didn't do these changes).

I'm loathe to go down a rabbit hole for a concept that seems to be a more natural fit for an entirely different system (FaceFaceConstraints)...

I asked for help in converting the FaceFaceConstraint (which now is deprecated and is renamed as Mortar) from mortar to DG a number of times. I also asked for suggestion on how to do it myself as you might have more compelling things to work on (email and this discussions here) but you know better than what the answer has been until today. If you can provide me with guidance in doing the conversion from mortart to DG or penalty I can do that (this is not a moose user change, is more of adding a core functions that requires deep knowledge of the framework). HoweverI have to say that conventional Cohesive elemnents does not differe so much in their limitation to the interface kernel where the neighboorhood is considered to be the one present in the undeformed configuration (i.e. calculation are performed always  between the same faces). These is why we decided to follow the route Wen implemented months ago. Also It was our only option at that time.

I hope this clarify your question and concerns.
Kind regards

Andrea
 

--

You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/29a2902a-3c90-4f81-90cf-008997fb5d70%40googlegroups.com.

[walkand...@gmail.com]

Hi andrea, 

a small issue in /materials/DisplacementJumpBasedCohesiveInterfaceMaterial.C

link: https://github.com/arovinelli/moose/blob/czm_DG_%2311456/modules/tensor_mechanics/src/materials/DisplacementJumpBasedCohesiveInterfaceMaterial.C

*******************************************************

 error: ‘ccoouutt’ is not a member of ‘std’

 std::ccoouutt << "==========================================================================="

*******************************************************

it should be cout right ?

best regards

在 2018年5月18日星期五 UTC+2下午5:47:37，Andrea Rovinelli写道：

[Andrea Rovinelli]

Yes you are right, I changed it to ccoouutt to make it pass moose tests.

However all those lines should be commented.

I’ll update the bar ache and squash it later in the morning

Sent from my iPhone

--
You received this message because you are subscribed to a topic in the Google Groups "moose-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/moose-users/hQkaOVbrm88/unsubscribe.
To unsubscribe from this group and all its topics, send an email to moose-users...@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/63c9afc9-1c8d-4746-9f5b-72cabb5c6bce%40googlegroups.com.

[walkand...@gmail.com]

Hi andrea, 

thank you sincerely for the contribution.

I read the code, and have a question about the following code in 'DisplacementJumpBasedCohesiveInterfaceKernel' for the interface kernel,

the code locate at the function for off-diag jacobian calculation:

Real DisplacementJumpBasedCohesiveInterfaceKernel::computeQpOffDiagJacobian

here:

***************************

if (jvar != _disp_1_var && jvar != _disp_2_var && jvar != _disp_1_neighbor_var &&

      jvar != _disp_2_neighbor_var)

  {

    mooseError("wrong variable requested");

  }

***************************

so if we have more dofs than ux,uy and uz, it seems this if-statement will always run the mooseError ?

for example ux,uy,uz with non-conserved order parameter eta, then the var for eta can't equals to any of these _disp_var right?

maybe it should be:

 if(jvar in offdiag_disp_var)

   do off-diag

 else if(jvar in other field)

   do something or return 0

 else

   mooseError('invalid dof var')

I'm not sure my understanding is correct or not ?

Thank you sincerely for the CZM work in MOOSE.

best regards

Yang

在 2018年5月23日星期三 UTC+2下午2:43:01，Andrea Rovinelli写道：

[Andrea Rovinelli]

Yang,
you are correct, but it is done on purpose. That specific interface kernel works for cohesive laws depending only on the displacement jump (u_slave-u_master).
If you need more variables first you need to couple them (you will need to use an additional kernel for each additional variable scalar component), then you will also need to update residuals and jacobian.
You can always start form that kernel and extend it to the problem you are facing. Note that you might also need to do some other work on the material and userobject. Just use the code provided as a base and keep adding complexity.
Hope this helps.

[walkand...@gmail.com]

Hi andrea, 

thanks for the reply.

Now I understand.

Another question relate to the development of new code.

I see you create InterfaceKernel and MeshSplit code under Tensor_Mechanics module,

I also try to add some code to the Tensor_Mechanics module, I just add the .h file to the related include dir,

and the .cpp file to its related src dir.

But when I try to run the 'make' command under Tensor_Mechanics folder, it seems it is not on the compile list.

So if I want to add something like you did for CZM, where should I modify the code and let it work?

Sorry for the poor experience in programming.

Best regards

在 2018年5月23日星期三 UTC+2下午7:38:32，Andrea Rovinelli写道：

[Andrea Rovinelli]

For this you may want to reference to the moose examples webpage http://mooseframework.org/wiki/MooseExamples/ .

Follow the examples and you will learn how to add moose objects.

The short answer is: put the .h file in include and .C (not cpp) file in src folder. 

Andrea

[Martin Lesueur]

Hi Andrea,

Thank you for your contribution. I am trying to do a similar InterfaceKernel and I had two questions about your implementation if you don't mind.

I was wondering how you are managing to keep only one variable of displacement between the two blocks. 

I understand that a necessary solution is to duplicate the nodes at the interface so that you can have a displacement discontinuity at the interface. This is why you created the CohesiveZoneMeshSplit. I do not understand however how this solution is sufficient; I do not understand how the code manages to distinguish between disp and disp_neighbor, since both variable and neighbor variable are disp in the InterfaceKernel block of your input file.

How does MOOSE know how to switch from Moose::Element to Moose::Neighbor?

When it will be clarified for me, I would like to know if you think that your CohesiveZoneMeshSplit is quite generic in the sense that it can be applied to all types of InterfaceKernel that are linking the same variable across the interface. Or is there some restrictions to using it?

I tried to use a mesh generated by your CohesiveZoneMeshSplit with my code and switch to only one variable in my input file and it didn't work. But it might be because I am using MatchValueBC that has access to a coupled_var and not a neighbor_var. I will try to implement MatchValueBC in an InterfaceKernel...

The second question can be extended to MOOSE developers as well. It regards the general way to implement InterfaceKernel.

I noticed that there exists two source files in MOOSE that are implementing the same equation:

InterfaceDiffusion in Moose/test

InterfaceDiffusionFluxMatch in the phase field module

The interesting thing is that the implementation differs between the two. We are comparing:

  Real r = 0.5 * (-_D[_qp] * _grad_u[_qp] * _normals[_qp] +
                  -_D_neighbor[_qp] * _grad_neighbor_value[_qp] * _normals[_qp]);

  switch (type)
  {
    case Moose::Element:

      r *= _test[_i][_qp];
      break;

    case Moose::Neighbor:
      r *= -_test_neighbor[_i][_qp];
      break;
  }

to

  Real res =
      _D * _grad_u[_qp] * _normals[_qp] - _D_neighbor * _grad_neighbor_value[_qp] * _normals[_qp];

  switch (type)
  {
    case Moose::Element:

      return res * _test[_i][_qp];

    case Moose::Neighbor:
      return res * _test_neighbor[_i][_qp];
  }

I personally have chosen the latter while Andrea chose to implement a modified version of the former (because Andrea is using the same variable across the interface?) 

Is there any difference to implement one or the other?

Thank you in advance for taking the time to answer my questions.

Cheers,

Martin

[Martin Lesueur]

Update before everyone wakes up in America, changing the MatchValueBC to an InterfaceKernel did not help.

Also is it possible for you Andrea to explain or give a reference for the way you computed the local traction in computeTractionLocal from CohesiveLaw_3DC?

It does not match with the mechanical traction that I computed.

Cheers

[Alexander Lindsay]

This one immediately above gives wrong results. Try running it on the `coupled_value_coupled_flux.i` input file, and you will see that the intersection point is at .8. The analytically correct answer is at .67 (2/3).

Alex

 

--

You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/f0b663ac-d906-4c5d-b230-fde8a7ebaaa4%40googlegroups.com.

[Andrea Rovinelli]

Dear Martin,
first let me say I'm still working on this so don't take the code as perfect at the moment.

I was wondering how you are managing to keep only one variable of displacement between the two blocks.

I understand that a necessary solution is to duplicate the nodes at the interface so that you can have a displacement discontinuity at the interface. This is why you created the CohesiveZoneMeshSplit. I do not understand however how this solution is sufficient; I do not understand how the code manages to distinguish between disp and disp_neighbor, since both variable and neighbor variable are disp in the InterfaceKernel block of your input file.

 The displacement field is described by one variable for each coordinate (disp_x, disp_y, disp_z) at each node. If you split the mesh by adding additional nodes, you will have a different value of disp_x, disp_y, disp_z on each side of the interface, however they will still be represented by the same variable_name. For a 1D case you start from

0----0----0 with node ids
1     2     3

and end up with
  
0----0 0----0 with node ids
1     2 3     4

now the variable for node 2 and 3 is the same (i.e. disp_x) but their value can be different and you can compute a jump like (disp_x_2 - disp_x_3)

How does MOOSE know how to switch from Moose::Element to Moose::Neighbor?

You don't really need to switch between them in the case of cohesive zone model, because the residual is going to be equal and opposite on both sides of the interface because of the continuity of traction derived from he equilibrium equation. The important thing is to choose a reference side and compute the displacement jump accordingly (i.e.  if reference side is the one with node 2 then D=disp_x_3 - disp_x_2 and viceversa). With this in hand you compute the residual accordingly to you traction separation law and set it positive on the reference and negative on the other side.

For the jacobian is the same, as your residual will always be equal and opposite by construction (however you have to take care of 4 cases).

Note that I need a spitted mesh because I want to see separations of blocks. You could have achieved the same results by using two displacement fields one for each side with monolithic mesh.
However this would results in a very complicated input file if you have many blocks with many interfaces.

When it will be clarified for me, I would like to know if you think that your CohesiveZoneMeshSplit is quite generic in the sense that it can be applied to all types of InterfaceKernel that are linking the same variable across the interface. Or is there some restrictions to using it?

I tried to use a mesh generated by your CohesiveZoneMeshSplit with my code and switch to only one variable in my input file and it didn't work. But it might be because I am using MatchValueBC that has access to a coupled_var and not a neighbor_var. I will try to implement MatchValueBC in an InterfaceKernel...

The CohesiveZoneMeshSplit does what I illustrated before but is not a mesh modifier. It needs an exodus mesh file with blocks already defined in it. i don't think you can uset on amoose generated mesh.
It is quite generic, the only limitation is the number of blocks and boundaries that  it can handle, but is the same that MOOSE can handle.
Ideally it could be used  in any situation working on an interfaceKernel with advantage being that you need only one variable on the entire domain.

The second question can be extended to MOOSE developers as well. It regards the general way to implement InterfaceKernel.

I noticed that there exists two source files in MOOSE that are implementing the same equation:

InterfaceDiffusion in Moose/test

InterfaceDiffusionFluxMatch in the phase field module

The interesting thing is that the implementation differs between the two. We are comparing:

  Real r = 0.5 * (-_D[_qp] * _grad_u[_qp] * _normals[_qp] +
                  -_D_neighbor[_qp] * _grad_neighbor_value[_qp] * _normals[_qp]);

  switch (type)
  {
    case Moose::Element:
      r *= _test[_i][_qp];
      break;

    case Moose::Neighbor:
      r *= -_test_neighbor[_i][_qp];
      break;
  }

to

  Real res =
      _D * _grad_u[_qp] * _normals[_qp] - _D_neighbor * _grad_neighbor_value[_qp] * _normals[_qp];

  switch (type)
  {
    case Moose::Element:
      return res * _test[_i][_qp];

    case Moose::Neighbor:
      return res * _test_neighbor[_i][_qp];
  }

I personally have chosen the latter while Andrea chose to implement a modified version of the former (because Andrea is using the same variable across the interface?) 

Is there any difference to implement one or the other?

 The second seems to be wrong as the residual is the same on both sides (there should be a minus) somewhere

Hope this clarify your questions and helps you going further.

Best

Andrea

[Martin Lesueur]

Thanks for your answers Andrea and Alex,

I think that first of all I should introduce my InterfaceKernel.

It is solving for the continuity of each component of the traction. The traction is defined in mechanics by stress*normal_vector.

So my implementation is:

    // continuity of each component of traction: sigma_ij . n_i

    RankTwoTensor stress_diff = _stress_master[_qp] - _stress_slave[_qp];

    RealVectorValue res = stress_diff * _normals[_qp];

    switch (type)

    {

      case Moose::Element:

        return res(_component) * _test[_i][_qp];

      case Moose::Neighbor:

        return res(_component) * _test_neighbor[_i][_qp];

This InterfaceKernel has to be coupled with MatchedValueJumpBC, deriving from MatchedValueBC where the residual is var - coupled_var - delta_disp.

This implementation is working, I have one benchmark for example where I have no delta_disp and I validate the solution with interface to the continuous problem which should be the same. With the interface, the traction is continuous and is equal to the value of the reference. The two files are compression and compression_interface.

After reading your posts, I tried the other implementation to see if my tests are still passing and they are!

This new implementation looks like that:

    // continuity of each component of traction: sigma_ij . n_i

    RankTwoTensor stress_diff = - _stress_master[_qp] - _stress_slave[_qp];

    RealVectorValue res = stress_diff * _normals[_qp];

    switch (type)

    {

      case Moose::Element:

        return res(_component) * _test[_i][_qp];

      case Moose::Neighbor:

        return - res(_component) * _test_neighbor[_i][_qp];

Note that if I multiply the residual by 0.5 (like in InterfaceDiffusion) it does not work anymore.

I am now confused since you told me that the first implementation (InterfaceDiffusionFluxMatch) is supposed to be wrong. However I know it is giving me the right results. And I know now as well that the implementation you are vouching for (InterfaceDiffusion) is also right.

However I do not really understand the implementation of InterfaceDiffusion. Is it possible to clarify how does Moose process this residual, which contribution is given to each residual and how does it end up reconciling to the final physical equation which the flux continuity? That would be really helpful, thanks.

Coming back on the mesh splitting, I think it's a great addition to MOOSE framework Andrea, so I tried to make it work for the Moose test of the InterfaceKernel, coupled_value_coupled_flux. I created a generated mesh from Moose that I then opened it with your CohesiveZoneMeshSplit that created the mesh for my input file: mesh_duplicate-nodes (see attached). I then ran it with a modified version of the test file (see attached), that uses now only one variable but the simulation did not give me the right results...

Is it possible that you have a look Andrea to see if it is possible to make it work? I would love to see this functionality in the MOOSE framework because I think it would be useful to a few people using InterfaceKernels.

I think it would also alleviate the work on your current pull request where you could take the mesh splitting out in a separate pull request to add in the MOOSE framework (and not the TensorMechanics). Let me know your thoughts about this.

Cheers,

Have a good week-end,

Martin

[Andrea Rovinelli]

Dear Martin,
First thanks for having fun with the CohesiveZoneMeshSplit, it's always nice when someone takes at good use your work.

Thanks for your answers Andrea and Alex,

I think that first of all I should introduce my InterfaceKernel.

It is solving for the continuity of each component of the traction. The traction is defined in mechanics by stress*normal_vector.

So my implementation is:

    // continuity of each component of traction: sigma_ij . n_i

    RankTwoTensor stress_diff = _stress_master[_qp] - _stress_slave[_qp];

    RealVectorValue res = stress_diff * _normals[_qp];

    switch (type)

    {

      case Moose::Element:

        return res(_component) * _test[_i][_qp];

      case Moose::Neighbor:

        return res(_component) * _test_neighbor[_i][_qp];

This InterfaceKernel has to be coupled with MatchedValueJumpBC, deriving from MatchedValueBC where the residual is var - coupled_var - delta_disp.

This implementation is working, I have one benchmark for example where I have no delta_disp and I validate the solution with interface to the continuous problem which should be the same. With the interface, the traction is continuous and is equal to the value of the reference. The two files are compression and compression_interface.

After reading your posts, I tried the other implementation to see if my tests are still passing and they are!

This new implementation looks like that:

Good to hear that :)
 

    // continuity of each component of traction: sigma_ij . n_i

    RankTwoTensor stress_diff = - _stress_master[_qp] - _stress_slave[_qp];

    RealVectorValue res = stress_diff * _normals[_qp];

    switch (type)

    {

      case Moose::Element:

        return res(_component) * _test[_i][_qp];

      case Moose::Neighbor:

        return - res(_component) * _test_neighbor[_i][_qp];

Note that if I multiply the residual by 0.5 (like in InterfaceDiffusion) it does not work anymore.

I am now confused since you told me that the first implementation (InterfaceDiffusionFluxMatch) is supposed to be wrong. However I know it is giving me the right results. And I know now as well that the implementation you are vouching for (InterfaceDiffusion) is also right.

However I do not really understand the implementation of InterfaceDiffusion. Is it possible to clarify how does Moose process this residual, which contribution is given to each residual and how does it end up reconciling to the final physical equation which the flux continuity? That would be really helpful, thanks.

If you really want understand this you may wnat to look at the InterfaceKernel.C and .h. Alex is the best person for answering this question.My answer would be
The Implementation really depends upon your problem. The 0.5 is intended to compute an average value and it might not the proper way of imposing your residual.
As long as your math is consistent with your reference frame you should be fine. Just make sure that your satisfying the equilibrium equation (i.e summing traction on both side of the interface should give you zero). The proper sign of the residual can be achieved by using either the one normal and changing the sign of the stresses on one side, or using the both normals and not flipping the stress tensor. This is up to you and is really problem dependent (i.e. if you have large deformation you may want use both normals and not flipping the stress).
I have to say that imposing the residual on the traction itself is somehow unusual for a cohesive model, usually I impose the residual on the displacements, because their redisual is the nodal force (i.e. the traction). This is just a thought anyway.
 

Coming back on the mesh splitting, I think it's a great addition to MOOSE framework Andrea, so I tried to make it work for the Moose test of the InterfaceKernel, coupled_value_coupled_flux. I created a generated mesh from Moose that I then opened it with your CohesiveZoneMeshSplit that created the mesh for my input file: mesh_duplicate-nodes (see attached). I then ran it with a modified version of the test file (see attached), that uses now only one variable but the simulation did not give me the right results...

I'm not sure what you mean by is not giving you the right results? Could you please tell us exactly which test are you trying to match? There is no coupled_value flux with the same mesh geometry you posted. Also looking at your input file there is one mistake.
You need to call the mesh split function explicitly on a monolithic mesh and there is a reason for doing this. If you start with a splitted mesh, the internal mapping of element neighbor is lost.

[Mesh]
 type = FileMesh
 file = mesh_duplicate-nodes.e
[]

Instead if you do something like below (note that I'm missing some parameter here)the mesh will be splitted but the element neighbor map will be preserved, does making the interfacekernel work properly.

[Mesh]
 type = CohesiveZoneMeshSplit
 file = mesh2split.e
[]

Try to do this and then rerun your example and let us know if results are as expected.

 

Is it possible that you have a look Andrea to see if it is possible to make it work? I would love to see this functionality in the MOOSE framework because I think it would be useful to a few people using InterfaceKernels.

I think it would also alleviate the work on your current pull request where you could take the mesh splitting out in a separate pull request to add in the MOOSE framework (and not the TensorMechanics). Let me know your thoughts about this.

I didn't thought about using the meshSplit as general tool to define only one variable. Let's see after we make this test work.

Have a nice weekend

Andrea
 

[Martin Lesueur]

Hey guys,

Thanks to your correction Andrea and the realisation that all other items relating to the interface such as MatchedValueBC and FluxBC need to be implemented as InterfaceKernels with your mesh splitting, I have managed to make the test work! The two input files (to be run with MooseTest using Alex' pull request and the mesh splitting of Andrea, CohesiveZoneMeshSplitBase and CohesiveZoneMeshSplit) are attached and give the same results.

I have also successfully switched to one variable in my work and still validated my previous benchmarks!

As I said though some items could not be used as before:

MatchedValueBC needs to become PenaltyInterfaceDiffusion (from Alex' new pull request)

DiffusionFluxBC needs to become something like InterfaceDiffusionBoundaryTerm from Phase Field

With this working, I would really like to push this feature in MOOSE framework, ideally as a meshmodifier that takes in input a sideset and duplicates its node to be used as a interface.

Do you think that could be done, Alex, Andrea?

I hope we can make it work, I'm really happy in any case for my research ^^ Thanks Andrea!

[Andrea Rovinelli]

Dear Martin,

I'm glad to hear you have been able to use my mesh splitting code to reduce the number of variables.

Implementing the mesh split as mesh modifier it's possible and it's something that I should have done from the beginning.

However, implementing the mesh splitting by sidesets rather than by blocks is not trivial, because there are cases in which it will fail. For  instance assume a 3x3 quad4 mesh like in the image below, in which you want to split internal sides of the center element only for faces aligned with the y normal. In this case I can't think of a way to split the mesh. There are geometric restrictions that I need to figure out (i.e. never duplicate nodes at the edge of the sidesets if they are in the bulk of the mesh).

I will add this to my to do list.

As for now I can modify my code to work as mesh modifier, submit a PR and make it available to the community quickly.

If you have any suggestions, you are more than welcome to post them in the here Add a mesh modfier to split a mesh by duplicating nodes.

I will open a PR for the mesh split in the next few days as soon as the code is translated from mesh to mesh modifier (deadline are coming so don't be surprised if the actual implementation and deployment to MOOSE is delayed).

Cheers

Andrea

[Andrea Rovinelli]

Martin,

I was looking to your input file, why do you need

[MeshModifiers]

  [./break_boundary]

    type = BreakBoundaryOnSubdomain

  [../]

[]

The CohesiveZoneMeshSplit should already take care of everything.

Andrea

[Martin Lesueur]

The Meshmodifier was just put in order to compare two input files that are as identical as possible. But we should remove it indeed because it is showing one advantage of the mesh splitting.

Great to hear that you are working on it. Having the input as the two blocks that separate the interface instead of the sideset is not a problem at all, I think that all of the interfacekernels users have two different blocks separating their interface in any case.

Cheers,

Martin

[Martin Lesueur]

Actually, I think that MatchedValueBC is not the same as PenaltyInterfaceDiffusion. The first one is enforced while the second one is added to the residual. At least for me, switching from one to the other did not give me the same results, even when playing with the penalty value.

However we can not use MatchedValueBC with the mesh splitting and this is probably because we access coupled_car instead of neighbor_var. I have also tried to write the BC in an InterfaceKernel (like an "InterfaceBC") without multiplying the condition by the test function but that did not work either.

Maybe the solution is to actually create a real class InterfaceBC?

This BC is used in the interfacekernel tests of moose so it would be good to have a compatibility with this new functionality.

What do you think?

[Alexander Lindsay]

Sounds like continuous creep of the InterfaceKernel system into the Constraint system realm...

--

You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/a8a0334d-b39a-40bc-8337-0bbe4569d0c6%40googlegroups.com.

[walkand...@gmail.com]

Hi Andrea,

  sorry to disturb you.

  I tried your mesh split code for several tests.

  I found that, sometime this mesh split code can't split the interface completely.

  I tried it on a plane domain, which looks like:

only two blocks, interface is the line between two plane.

I used peacock to check the split mesh , and it looks like:

for the nodesets:

for the boundaries:

it looks like some interface node information is lost.

I applied a flux at the right side of the plane for a diffusion problem, results show that some part of the interface is not splitted, which means

when I apply the continue condition on the interface with an interface kernel(just simple penalty method to let c on the interface to be continued), 

the result is not continue on this part, which looks like:

So I'm wondering how can I solve this problem?

Thank you 

Best regards

Yang

在 2018年7月1日星期日 UTC+2下午6:36:47，Andrea Rovinelli写道：

[Andrea Rovinelli]

Dear Yang, thanks for letting me know.

First I have one question, are you using the latest version of my code (the one working as a mesh modifier) or the previous one (the one working as a Meshreader and modifier)?

In either case, could you please send me or post the mesh file that is causing this problem.

I would like to dig a bit more on this issue.

Thanks.

Andrea

[walkand...@gmail.com]

Hi Andrea,

  thanks for the reply.

  I think I'm still using the old version. I will try the new one from you github repository.

  

  The attachment is the geo files, one is the plane domain, another one for a circle with several domains .

  Thank you sincere

  Best regards

  Yang

在 2018年7月7日星期六 UTC+2下午4:21:58，Andrea Rovinelli写道：

[Andrea Rovinelli]

Dear Yang,

the ones you attached are geometry files (.geo), could you post the .msh or whatever format you are exporting to and use in your MOOSE input file?

Thanks

Andrea

[walkand...@gmail.com]

Hi Andrea,

 sorry for the mistake.

 The attachment is the .msh file.

  By the way, where can I find these CohesiveZoneMeshSplit.cpp CohesiveZoneMeshSplit.h things?

  I git cloned your czm_reworking branch, it seems they are not there?

  I remember in the old version, for both the include and src path inside the tensor_mechanics module, there is a folder named mesh

  contains the .h and .cpp file.

  Best regards 

  Yang

在 2018年7月7日星期六 UTC+2下午4:36:01，Andrea Rovinelli写道：

[Andrea Rovinelli]

Yang, Mesh split is now called BreakMeshByBlock and you can find the relevant files into moose/framework/src/meshmodifiers and moose/framework/include/meshmodifiers.

Also, please make sure to run your mesh as replicated, add the following line into your mesh block in the input file

  parallel_type = replicated

Right now this tool does not support distributed mesh (it should error out if you don set it correctly).

Thanks

[Martin Lesueur]

Hi Andrea,

I confirm that the meshes from gmsh were not splitted properly with your last version of the meshsplitting.

I believe it will be different when using it as a meshmodifier but we shall see with Yang's new results.

[walkand...@gmail.com]

Hi Andrea and Martin,

  I tried the new 'BreakMeshByBlock' things, for the mesh part, my input file looks like:

[Mesh]

  type = FileMesh

  file = plane2.msh

  parallel_type = REPLICATED

[]

[MeshModifiers]

  [./meshsplit]

    type = BreakMeshByBlock

    boundaries = inner

    interface_id = 5001

  [../]

[]

these incomplete things seem still there.

Here is my result:

I also tried to use different mesh(by change the 'dx' parameter in my geo file, which can gives me the mesh from coarse case to fine case)

so for dx=0.4, only one node is lost

for dx=0.3, two nodes are lost:

for dx=0.2, three line segment(not node) are lost:

for dx=0.1, five line segments are lost:

for dx=0.05, seven line segments are lost:

Best regards

Yang

在 2018年7月8日星期日 UTC+2上午3:58:47，Martin Lesueur写道：

[Andrea Rovinelli]

Dear Yang and Martin,

thanks for the patience.

I've fixed the code you my want to try to pull the new version from here

https://github.com/arovinelli/moose/tree/MeshSplit

I tried with both the meshes you shared and it works fine now.

Thanks Andrea

[walkand...@gmail.com]

Thank you so much Andrea,

yes, it works now.

Thank you .

Best regards

Yang

在 2018年7月9日星期一 UTC+2下午8:21:01，Andrea Rovinelli写道：

[Andrea Rovinelli]

Great :)

[Martin Lesueur]

Hi Andrea,

Congrats on the pull request!

I'm sad though because I can't use the mesh splitting just yet. It is not compatible with my interface implementation.

I want to enforce traction continuity and displacement continuity.

To do so I was using MatchValueBC at the interface which links disp0 to disp1. For the traction, I implemented an interfacekernel that makes the traction of stress0 and stress1 equal.

This was working well.

If I switch to one variable with the mesh split, I still have the continuity of stress because of the compatibility of the interfacekernels with your mesh split. But I cannot reproduce the displacement discontinuity because the MatchValueBC takes a coupled_variable and does not recognize a neighbor_var.

Any idea how I could make it work?

For your implementation Andrea, how did you manage to do it?

Cheers,

Martin

[Martin Lesueur]

To ask a more specific question, is it possible to access neighbor_var for a BC?

[Alexander Lindsay]

For an IntegratedBC you cannot use the values or gradients from a variable that exists only on the neighboring subdomain

--

You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/bc1e1c82-f11e-4851-bd62-4668dbfdd20b%40googlegroups.com.

[Andrea Rovinelli]

Dear Martin,
I believe MatchedValueBC is not the appropriate object for your problem (it impose the value of two variables to match on a boundary, e.g. on side).
Why don't you impose the continuity of traction with another interface kernel? However, I don't get why you need to impose the continuity of traction is it should be automatically granted by the enforcing equilibrium.

To answer your question. In my model the traction at the interface is a function of the displacement jump. To make it simple T=K (u+ - u-).  The residual of the displacements are forces. on the interface this forces are tractions. Therefore my interface kernel impliclty impose a residual equal to the tractions, which are calculated depending on the displacement jump.

Continuity of displacement means  u+ - u- = 0. So, I guess your interface kernel does something like R = (u+ - u-) K, where K is a penalty tern, for each component of the displacement vector.
However you can also interpret K as the stiffness of your interface (your case fall back to my implementation). In my mind you don't need to enforce traction continuity as you are already doing it implicitly.
To check this is what is happening behind the scene, declare 2/3 nodal aux variables (F_x, F_y,F_z) and use "save_in='F_x F_y F_z' " in the Tensor mechanics block. This will save the value of the residual of displacement at each node (traction)
Compare the values of those values on both side of you interface and you will see that they match. Those are your traction.

Hope this help.

Andrea

[Martin Lesueur]

Alex,

The equivalent of my IntegratedBC is already working with the mesh split because it is implemented in my InterfaceKernel that sets the continuity of the traction.

But I am looking to enforce the MatchValueBC (which is a NodalBC) with the mesh split.

I have seen in this post that you are doing the same as me. So if you wanted to switch to one variable, how you would do as you say in your post:

The matched value bc strongly sets one variable, let's call it 'u', equal to the value of the other, let's call it 'v'

where v is now u of the neighbor node in the interface?

Andrea, I am sorry I have a hard time understanding your implementation... Do you have perhaps some literature that explains this more in detail?

Cheers,

Martin

I think it would also alleviate the work on your current <a href="https://github.com/idaholab/moose/pull/11

[Martin Lesueur]

When you are desperate, it's not because you don't understand it that you won't try it XD

Andrea what you are saying is true, I ran my tests with PenaltyInterfaceDiffusion (which is the equivalent of MatchedValueBC as an InterfaceKernel) instead of MatchedValueBC and I deactivated my InterfaceKernel for the traction continuity and it converged to the right solution! There was a small exodiff but i'm happy with it.

With this implementation I can finally switch to one variable using your mesh split!

So now, I REALLY want to understand the details of your derivation. As I said in my previous post, if you have any literature that I could refer to, it would be really useful...

Don't tell my supervisors but I'm still happy it works even though I don't understand it ;)

[Andrea Rovinelli]

Martin,
please find attached a paper about discontinuous galerkin method for cohesive zone modeling.

On Friday, July 20, 2018 at 1:36:19 AM UTC-5, Martin Lesueur wrote:

When you are desperate, it's not because you don't understand it that you won't try it XD

:)

 

Andrea what you are saying is true, I ran my tests with PenaltyInterfaceDiffusion (which is the equivalent of MatchedValueBC as an InterfaceKernel) instead of MatchedValueBC and I deactivated my InterfaceKernel for the traction continuity and it converged to the right solution! There was a small exodiff but i'm happy with it.

The small difference in the solution is unavoidable using this method, you are not enforcing exactly the same displacement, you are using a penalty to enforce it (you might want to calibrate your penalty value somehow).

With this implementation I can finally switch to one variable using your mesh split!

So now, I REALLY want to understand the details of your derivation. As I said in my previous post, if you have any literature that I could refer to, it would be really useful...

What you are doing is described in the paper as  intrinsic TSL (Fig 1 (a)).
The value you assign to the PenaltyInterfaceKerrnel is K in Figure 1(a), the only difference is that you assume an interface with linear behavior (you just use the linear portion).
 

Don't tell my supervisors but I'm still happy it works even though I don't understand it ;)

I won't ;)

Hope this help.

Andrea

 

</block

[walkand...@gmail.com]

Hi Andrea,

Sorry to disturb you, I find a new issue for the BreakMeshByBlock function.

Maybe it is just a small bug...since in most cases, it works very well.

I tried to apply a simple shear load to a rectangle domain(just simple linear elastic problem, without any cohesive law or interface kernels).

First case: apply shear loading(let disp_y=0.01*t and -0.01*t on the top-left and bottom-left edge respectively), at the beginning it looks like:

after several steps, it looks like:

For case two: I apply the same boundary condition but on the right side

at the beginning it looks like:

after several steps, it looks like:

I have no idea how this happens, it seems inside the interface, everything works fine, but the problem happens 

for the two nodes at the end side?

The attachment is my input file and mesh file.

If you have time could help me to figure out this issue?

Question 2:

Is it possible to store an interface material property and history material property?

I try to implement the PPR cohesive model(https://doi.org/10.1016/j.engfracmech.2012.02.007), the model presented by Park where the maximum history displacement jump is required for the calculation.

I try to define and also get material property old stuff, but it seems not work: https://groups.google.com/forum/#!searchin/moose-users/walkandthinker|sort:date/moose-users/MOGZPwvPgxY/Vi15uMFZCQAJ

Best regards

Yang

在 2018年2月1日星期四 UTC+1下午11:53:07，Andrea Rovinelli写道：

Dear all,
I need to setup a cohesive zone model to study intra-granular cracking. In this model cracking can occur only at grain boundaries and to implement this behavior I was thinking to use cohesive elements. I would have XFEM because I don't need an arbitrary crack path, but just grain boundary cracking.
First I'm not sure if MOOSE can handle cohesive elements, i tried to look into libmesh supported element types (https://libmesh.github.io/doxygen/classlibMesh_1_1Elem.html) and I don't see specific elements allowed to have zero thickness. In fact when I try a simple 2D case in which I add some zero thickness QUAD4 elements moose complains saying that the element has zero volume.

I know someone from the MOOSE development team is working on this feature and that at some point it will be released, but we need to work on this project.

Any hint and help would be really appreciated.

Thanks in advance


Andrea (ANL)

[Wen Jiang]

Hi Yang,

We recently found several issues with the BreakMeshByBlock. It seems one of your issue is related to that. I opened a ticket here https://github.com/idaholab/moose/issues/12033 and I will be working on it. Thanks for pointing out the issue. 

Wen