Feature nonlinear #3

kianoosh1989 · 2025-04-08T16:23:44Z

Dear Dr. Najian,

The following PR contains the paper's nonlinear 2D mechanic example. Additionally, some modifications were necessary to account for nonlinear Neo-Hookean material.

Kind regards,
Kianoosh

…odifications were implemented to account for nonlinearity.

RezaNajian · 2025-05-20T03:58:51Z

@kianoosh1989 Thanks for the PR. I merged the main to update the branch, however two tests are failing, please check and correct. It is due to the changes in the nonlinear solver.

RezaNajian · 2025-05-20T04:00:39Z

fol/loss_functions/mechanical_neohooke.py

Please make the loss function applicable to other elements types, like tri,tetra and so on. also update the integration test file to cover all elements.

RezaNajian · 2025-05-20T04:03:00Z

fol/tools/fem_utilities.py

        raise NotImplementedError("This method should be implemented by subclasses.")

    def fourth_order_identity_tensor(self, dim=3):
        """


Do not use for loops and jit the functions to improve the performance.

RezaNajian · 2025-05-20T04:03:57Z

fol/tools/fem_utilities.py

        """
-        C = jnp.zeros((dim, dim, dim, dim))
-        C = 0.5*(jnp.einsum('ik,jl->ijkl',A,B) + jnp.einsum('il,jk->ijkl',A,B))
-        return C


the same here do not use explicit for loop and reformulate them with other functionalities or hard code the set

RezaNajian · 2025-05-20T04:05:09Z

fol/tools/usefull_functions.py

    L = 1
    N = int((data.reshape(-1, 1).shape[0])**0.5)
    nu = loss_settings["poisson_ratio"]
    e = loss_settings["young_modulus"]


Please avoid making specific functions here move this to the example or where it is used. We need to clean this file.

RezaNajian

@kianoosh1989 I reviewed you PR and there are some points before we merge. Please apply and request again. Thanks for the good job.

RezaNajian

@kianoosh1989 Thankls a lot for the great work, for the sdake of maintainilty and cleanes i made some comments and suggestions let address them as much as possible before mergeing

RezaNajian · 2025-05-25T07:51:49Z

.gitignore

 *.zarr
 py-env*
+.venv
+*.txt


*.txt is too general to be in gitignore, please remove it from the list

RezaNajian · 2025-05-25T07:53:38Z

examples/mechanical_box/neo_hooke_mechanical.py

 import sys
 import os
-
+sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..','..')))


why do we need this add, it works without as well, if it is the case remove it

RezaNajian · 2025-05-25T07:55:03Z

examples/paper_iFOL/static_box/fourier_mechanical_implicit_pr.py

 # import necessaries 
 import sys
 import os
+sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..','..','..')))


the same goes here, please remove this append since it is redundant

RezaNajian · 2025-05-25T07:56:08Z

examples/paper_iFOL/static_box/mechanical3D_usefull_functions.py

+from fol.mesh_input_output.mesh import Mesh
+import copy
+
+def box_mesh(Nx, Ny, Nz, Lx, Ly, Lz, case_dir):


this function is general and can stay in the usefull function of the main. otherwise iot gets repeated in every box problem, please move it back to where it was

RezaNajian · 2025-05-25T07:57:53Z

examples/paper_iFOL/static_box/mechanical3D_usefull_functions.py

+
+    return meshio_obj
+
+def create_3D_box_mesh(Nx,Ny,Nz,Lx,Ly,Lz,case_dir):


maybe it was my mistake and I got it wrong but as I reciommended, please keep general functions like sample generatin and mesh genberation functions in useful functions of fol itself. Sorry for the confusion.

RezaNajian · 2025-05-25T08:12:43Z

fol/solvers/fe_nonlinear_residual_based_solver.py

    @print_with_timestamp_and_execution_time
-    def Solve(self,current_control_vars,current_dofs):
+    def Solve(self,current_control_vars,current_dofs_np:np.array):
+        current_dofs = jnp.array(current_dofs_np)


this line is redundant since you mentioned the type of input is numpy array, please update here

The reason for this casting from numpy to jnp.numpy is that I wanted to keep the main script appearance consistent for linear and non-linear cases. In other words, we give a zero numpy array to the FE solver in every main script as the initial guess, while for the nonlinear case, this array will be modified at each load step, and this modification happens in methods that take jax.numpy array types and return jax.numpy array types.

RezaNajian · 2025-05-25T08:15:45Z

fol/solvers/fe_nonlinear_residual_based_solver.py

                else:
                    fol_info(f"iteration:{i+1},delta_norm:{delta_norm},residuals_norm:{res_norm}")
-            current_dofs = applied_BC_dofs
+            current_dofs = current_dofs.at[self.fe_loss_function.non_dirichlet_indices].set(applied_BC_dofs[self.fe_loss_function.non_dirichlet_indices])


i am doubtfull about this line, could you please explain me wha tis happening here, teh computed deltas and dofs always satisfy the dirichlet BCs, and if they get added or subtracted they still satisfy or i do get the point here.

I just wanted to ensure that in every load step, after the solution has been computed for that step, only non-Dirichlet boundary values are replaced. Boundary conditions are definitely modified at each load step by the "ApplyDirichletBCOnDofVector" method, and only non-Dirichlet values are expected to be recalculated.

RezaNajian · 2025-05-25T08:18:54Z

fol/tools/fem_utilities.py

        xsie_vol = (k/4)*(J**2 - 2*jnp.log(J) -1)
        I1_bar = (J**(-2/3))*jnp.trace(C)
        xsie_iso = 0.5*mu*(I1_bar - 3)
-        loss_positive_bias = 100     # To prevent loss to become a negative number


what is going on here ?! i would not hard code sth like 100 for teh sake iof loss function sign. Could we come up with a cleaner solution ? a explanantion would be highly appreciated. here in this file we are doing FEM and should not care about training. right !

This part belongs to an old version that has been merged into the main branch and no longer exists in the feature_nonlinear. Please use the branch feature_nonlinear as the reference for any neo-Hookean material implementation.

tests/integration/test_mechanical_3D_hexa_poly.py

RezaNajian · 2025-05-25T08:23:30Z

tests/integration/test_neo_hooke_mechanical_3D_hexa_poly.py

the same comment for neo hook element, lets merge tests for diffrent element geometry types into a single file

RezaNajian · 2025-05-26T03:25:11Z

@kianoosh1989 I think you did not push the latest changes, please check

kianoosh1989 · 2025-05-26T09:24:05Z

All comments were applied in the latest push. Please kindly check the branch.

RezaNajian

@kianoosh1989 thanks for the updates, i have made some critical points regarding for loops in jitted functions. Please check. Also I noticed unit tests for neo hook element do not exist anymore, we need unit test for new core functionalities like element. Please bring them bacvk and we do not need integrations tests for all element types, only one is enough. I got confused. Please also add your examples to the example runner and make sure they get ran. Thanks a lot

RezaNajian · 2025-05-27T03:59:41Z

fol/loss_functions/mechanical_neohooke.py

-        gp_geo_stiffness = gp_geo_stiffness.at[4:6,2:4].set(jnp.eye(2) * (DN_DX_T[:,2].T @ (S_mat @ DN_DX_T[:,1])))
-        gp_geo_stiffness = gp_geo_stiffness.at[4:6,4:6].set(jnp.eye(2) * (DN_DX_T[:,2].T @ (S_mat @ DN_DX_T[:,2])))
-        gp_geo_stiffness = gp_geo_stiffness.at[4:6,6:8].set(jnp.eye(2) * (DN_DX_T[:,2].T @ (S_mat @ DN_DX_T[:,3])))
+        # Rearrange blocks into full (8, 8) matrix


These for loops degradate the performance since this function is going to be used in a computational intensiveenviroment. Please reformalte it and get rid of for loops.

RezaNajian · 2025-05-27T04:01:04Z

fol/loss_functions/mechanical_neohooke.py

-        gp_geo_stiffness = gp_geo_stiffness.at[4:6,2:4].set(jnp.eye(2) * (DN_DX_T[:,2].T @ (S_mat @ DN_DX_T[:,1])))
-        gp_geo_stiffness = gp_geo_stiffness.at[4:6,4:6].set(jnp.eye(2) * (DN_DX_T[:,2].T @ (S_mat @ DN_DX_T[:,2])))
+        # Rearrange blocks into full (8, 8) matrix
+        gp_geo_stiffness = jnp.block([[blocks[i, j] for j in range(num_nodes)] for i in range(num_nodes)])


the same here, we should not use explicit for loops inside jitted function sicne they all get rolled out and become easily inefficient

RezaNajian · 2025-05-27T04:01:23Z

fol/loss_functions/mechanical_neohooke.py

-        gp_geo_stiffness = gp_geo_stiffness.at[9:12,6:9].set(jnp.eye(3) * (DN_DX_T[:,3].T @ (S_mat @ DN_DX_T[:,2])))
-        gp_geo_stiffness = gp_geo_stiffness.at[9:12,9:12].set(jnp.eye(3) * (DN_DX_T[:,3].T @ (S_mat @ DN_DX_T[:,3])))
+        # Rearrange blocks into full (12, 12) matrix
+        gp_geo_stiffness = jnp.block([[blocks[i, j] for j in range(num_nodes)] for i in range(num_nodes)])  


the same here

RezaNajian · 2025-05-27T04:01:44Z

fol/loss_functions/mechanical_neohooke.py

-        gp_geo_stiffness = gp_geo_stiffness.at[21:24,18:21].set(jnp.eye(3) * (DN_DX_T[:,7].T @ (S_mat @ DN_DX_T[:,6])))
-        gp_geo_stiffness = gp_geo_stiffness.at[21:24,21:24].set(jnp.eye(3) * (DN_DX_T[:,7].T @ (S_mat @ DN_DX_T[:,7])))
+        # Rearrange blocks into full (24, 24) matrix
+        gp_geo_stiffness = jnp.block([[blocks[i, j] for j in range(num_nodes)] for i in range(num_nodes)])  


RezaNajian · 2025-05-27T04:03:00Z

tests/integration/test_mechanical_3D_poly.py

+        U_FEM = np.array(linear_fe_solver.Solve(K_matrix,np.zeros(U_FOL.shape)))
+        l2_error = 100 * np.linalg.norm(U_FOL-U_FEM,ord=2)/ np.linalg.norm(U_FEM,ord=2)
+        self.assertLessEqual(l2_error, 10)
+        print(l2_error)


please avoid explicit prints, so remove it

RezaNajian · 2025-05-28T04:08:13Z

@kianoosh1989 Thanks a lot for the changes. I went through them and it is indeed clean and organized.
One thing is that test_mechanical_poly.py is failing on github VM, therfore it does not allow to merge it.
you can check the test log undeer check page of this PR and i put the error here, could you please fix it and then merge it when all tests are passed on github. the error

kianoosh1989 · 2025-05-29T09:58:28Z

@RezaNajian There was a misspelling in a directory name, which caused the error during testing. It has now been resolved.

RezaNajian

@kianoosh1989 Thanks a lot for the crrections and great work. Let's merge it.

* setup environments * paper example: 2D_mechanical_nonlinear was added. In addition, some modifications were implemented to account for nonlinearity. * cleaning * Comments were addressed * the tests integration dir was modified in accordance with the comments * cleaning in accordance with the comments * Geometric stiffness became concise * for-loop was omitted from class mechanical_neohooke * examples were added to examples_runner.py * some cleaning * adjusting rtol and atol for test assertions * test_mechanical_3D_poly.py was checked --------- Co-authored-by: RezaNajian <[email protected]>

kianoosh1989 and others added 4 commits February 25, 2025 11:09

setup environments

5702a60

paper example: 2D_mechanical_nonlinear was added. In addition, some m…

1e4dc50

…odifications were implemented to account for nonlinearity.

cleaning

d75be2f

Merge branch 'main' into feature_nonlinear

bf6e212

RezaNajian reviewed May 20, 2025

View reviewed changes

RezaNajian requested changes May 20, 2025

View reviewed changes

kianoosh1989 added 2 commits May 22, 2025 21:52

Comments were addressed

0637f15

branch was updated with main and unit test for neo_hookean was added

7ea27e8

RezaNajian reviewed May 25, 2025

View reviewed changes

kianoosh1989 added 2 commits May 25, 2025 19:58

the tests integration dir was modified in accordance with the comments

2308fab

cleaning in accordance with the comments

417d0f2

Geometric stiffness became concise

9190b97

RezaNajian requested changes May 27, 2025

View reviewed changes

kianoosh1989 added 4 commits May 27, 2025 18:34

for-loop was omitted from class mechanical_neohooke

e9df7f6

examples were added to examples_runner.py

48e4374

some cleaning

d2b4c16

adjusting rtol and atol for test assertions

15efb86

test_mechanical_3D_poly.py was checked

7346289

RezaNajian approved these changes May 29, 2025

View reviewed changes

RezaNajian merged commit e912ffc into main May 29, 2025
1 check passed


		return meshio_obj

		def create_3D_box_mesh(Nx,Ny,Nz,Lx,Ly,Lz,case_dir):

Feature nonlinear #3

Feature nonlinear #3

Uh oh!

Conversation

kianoosh1989 commented Apr 8, 2025

Uh oh!

RezaNajian commented May 20, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

RezaNajian left a comment

Choose a reason for hiding this comment

Uh oh!

RezaNajian left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

RezaNajian commented May 26, 2025

Uh oh!

kianoosh1989 commented May 26, 2025

Uh oh!

RezaNajian left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

RezaNajian commented May 28, 2025

Uh oh!

kianoosh1989 commented May 29, 2025

Uh oh!