diff --git a/README.md b/README.md
index c4491335..96c04db5 100644
--- a/README.md
+++ b/README.md
@@ -16,10 +16,10 @@ scaleupROM is mainly built upon [MFEM](https://mfem.org/) and [libROM](https://w
   - Poisson equation
   - Stokes flow
   - Steady Navier-Stokes flow
+  - Linear elasticity
 
 ## Features to be added
 
-- Linear elasticity
 - Nonlinear elasticity
 - Time-dependent Navier-Stokes flow
 
diff --git a/examples/linelast/CMakeLists.txt b/examples/linelast/CMakeLists.txt
index 3fdd7dd8..8a673b45 100644
--- a/examples/linelast/CMakeLists.txt
+++ b/examples/linelast/CMakeLists.txt
@@ -2,5 +2,13 @@
 #
 # SPDX-License-Identifier: MIT
 file(COPY linelast.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast)
+file(COPY linelast.simpleL.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast)
+file(COPY config/linelast.simpleL.h5 DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/config)
+file(COPY linelast.lattice.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast)
+file(COPY config/linelast.lattice.h5 DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/config)
+
 file(COPY meshes/beam-tri.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes)
 file(COPY meshes/beam-tet.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes)
+file(COPY meshes/joint2D.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes)
+file(COPY meshes/rod2D_H.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes)
+file(COPY meshes/rod2D_V.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes)
\ No newline at end of file
diff --git a/examples/linelast/config/linelast.lattice.h5 b/examples/linelast/config/linelast.lattice.h5
new file mode 100644
index 00000000..93da8e1e
Binary files /dev/null and b/examples/linelast/config/linelast.lattice.h5 differ
diff --git a/examples/linelast/config/linelast.simpleL.h5 b/examples/linelast/config/linelast.simpleL.h5
new file mode 100644
index 00000000..b3a2fcc1
Binary files /dev/null and b/examples/linelast/config/linelast.simpleL.h5 differ
diff --git a/examples/linelast/linelast.lattice.yml b/examples/linelast/linelast.lattice.yml
new file mode 100644
index 00000000..e0013c6c
--- /dev/null
+++ b/examples/linelast/linelast.lattice.yml
@@ -0,0 +1,75 @@
+main:
+  mode: single_run
+  use_rom: true
+  solver: linelast
+
+mesh:
+  type: component-wise
+  component-wise:
+    global_config: "config/linelast.lattice.h5"
+    components:
+      - name: "joint2D"
+        file: "meshes/joint2D.mesh"
+      - name: "rod2D_H"
+        file: "meshes/rod2D_H.mesh"
+      - name: "rod2D_V"
+        file: "meshes/rod2D_V.mesh"
+
+domain-decomposition:
+  type: interior_penalty
+
+discretization:
+  order: 1
+  full-discrete-galerkin: true
+
+visualization:
+  enabled: false
+  unified_paraview: true
+  file_path:
+    prefix: lattice_output
+
+parameterized_problem:
+  name: linelast_disp_lattice
+
+single_run:
+  linelast_disp_lattice:
+    rdisp_f: 1.0
+
+sample_generation:
+  maximum_number_of_snapshots: 400
+  component_sampling: false
+  file_path:
+    prefix: "linelast_lattice"
+  parameters:
+    - key: single_run/linelast_disp/rdisp_f
+      type: double
+      sample_size: 1
+      minimum: 1.0
+      maximum: 1.0
+
+basis:
+  prefix: "linelast_lattice"
+  #number_of_basis: 8
+  tags:
+    - name: "joint2D"
+      number_of_basis: 12
+    - name: "rod2D_H"
+      number_of_basis: 8
+    - name: "rod2D_V"
+      number_of_basis: 9
+  svd:
+    save_spectrum: true
+    update_right_sv: false
+  visualization:
+    enabled: false
+
+model_reduction:
+  rom_handler_type: mfem
+  subdomain_training: universal
+  linear_solver_type: direct
+  save_operator:
+    level: component
+    prefix: "linelast.rom_elem"
+  compare_solution:
+    enabled: true
+    
\ No newline at end of file
diff --git a/examples/linelast/linelast.simpleL.yml b/examples/linelast/linelast.simpleL.yml
new file mode 100644
index 00000000..3775600a
--- /dev/null
+++ b/examples/linelast/linelast.simpleL.yml
@@ -0,0 +1,68 @@
+main:
+  mode: single_run
+  use_rom: true
+  solver: linelast
+
+mesh:
+  type: component-wise
+  component-wise:
+    global_config: "config/linelast.simpleL.h5"
+    components:
+      - name: "joint2D"
+        file: "meshes/joint2D.mesh"
+      - name: "rod2D_H"
+        file: "meshes/rod2D_H.mesh"
+      - name: "rod2D_V"
+        file: "meshes/rod2D_V.mesh"
+
+domain-decomposition:
+  type: interior_penalty
+
+discretization:
+  order: 1
+  full-discrete-galerkin: true
+
+visualization:
+  enabled: false
+  unified_paraview: true
+  file_path:
+    prefix: lcantilever_output
+
+parameterized_problem:
+  name: linelast_disp_lcantilever
+
+single_run:
+  linelast_disp_lcantilever:
+    rdisp_f: 10.0
+
+sample_generation:
+  maximum_number_of_snapshots: 400
+  component_sampling: false
+  file_path:
+    prefix: "linelast_Lbeam"
+  parameters:
+    - key: single_run/linelast_disp/rdisp_f
+      type: double
+      sample_size: 6
+      minimum: 10.0
+      maximum: 10.0
+
+basis:
+  prefix: "linelast_Lbeam"
+  number_of_basis: 2
+  tags:
+    - name: "comp0"
+  svd:
+    save_spectrum: true
+    update_right_sv: false
+  visualization:
+    enabled: false
+
+model_reduction:
+  rom_handler_type: mfem
+  subdomain_training: universal
+  save_operator:
+    level: component
+    prefix: "linelast.rom_elem"
+  compare_solution:
+    enabled: true
diff --git a/examples/linelast/meshes/joint2D.mesh b/examples/linelast/meshes/joint2D.mesh
new file mode 100644
index 00000000..69eca923
--- /dev/null
+++ b/examples/linelast/meshes/joint2D.mesh
@@ -0,0 +1,129 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+2
+
+elements
+16
+1 3 0 1 6 5
+1 3 1 2 7 6
+1 3 2 3 8 7
+1 3 3 4 9 8
+1 3 5 6 11 10
+1 3 6 7 12 11
+1 3 7 8 13 12
+1 3 8 9 14 13
+1 3 10 11 16 15
+1 3 11 12 17 16
+1 3 12 13 18 17
+1 3 13 14 19 18
+1 3 15 16 21 20
+1 3 16 17 22 21
+1 3 17 18 23 22
+1 3 18 19 24 23
+
+boundary
+16
+1 1 0 1
+1 1 1 2
+1 1 2 3
+1 1 3 4
+3 1 21 20
+3 1 22 21
+3 1 23 22
+3 1 24 23
+4 1 5 0
+4 1 10 5
+4 1 15 10
+4 1 20 15
+2 1 4 9
+2 1 9 14
+2 1 14 19
+2 1 19 24
+
+vertices
+25
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_2D_P1
+VDim: 2
+Ordering: 1
+
+0 0
+0.25 0
+0 0.25
+0.25 0.25
+0.25 0
+0.5 0
+0.25 0.25
+0.5 0.25
+0.5 0
+0.75 0
+0.5 0.25
+0.75 0.25
+0.75 0
+1 0
+0.75 0.25
+1 0.25
+0 0.25
+0.25 0.25
+0 0.5
+0.25 0.5
+0.25 0.25
+0.5 0.25
+0.25 0.5
+0.5 0.5
+0.5 0.25
+0.75 0.25
+0.5 0.5
+0.75 0.5
+0.75 0.25
+1 0.25
+0.75 0.5
+1 0.5
+0 0.5
+0.25 0.5
+0 0.75
+0.25 0.75
+0.25 0.5
+0.5 0.5
+0.25 0.75
+0.5 0.75
+0.5 0.5
+0.75 0.5
+0.5 0.75
+0.75 0.75
+0.75 0.5
+1 0.5
+0.75 0.75
+1 0.75
+0 0.75
+0.25 0.75
+0 1
+0.25 1
+0.25 0.75
+0.5 0.75
+0.25 1
+0.5 1
+0.5 0.75
+0.75 0.75
+0.5 1
+0.75 1
+0.75 0.75
+1 0.75
+0.75 1
+1 1
diff --git a/examples/linelast/meshes/rod2D_H.mesh b/examples/linelast/meshes/rod2D_H.mesh
new file mode 100644
index 00000000..14c5390b
--- /dev/null
+++ b/examples/linelast/meshes/rod2D_H.mesh
@@ -0,0 +1,393 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+2
+
+elements
+64
+1 3 0 1 18 17
+1 3 1 2 19 18
+1 3 2 3 20 19
+1 3 3 4 21 20
+1 3 4 5 22 21
+1 3 5 6 23 22
+1 3 6 7 24 23
+1 3 7 8 25 24
+1 3 8 9 26 25
+1 3 9 10 27 26
+1 3 10 11 28 27
+1 3 11 12 29 28
+1 3 12 13 30 29
+1 3 13 14 31 30
+1 3 14 15 32 31
+1 3 15 16 33 32
+1 3 17 18 35 34
+1 3 18 19 36 35
+1 3 19 20 37 36
+1 3 20 21 38 37
+1 3 21 22 39 38
+1 3 22 23 40 39
+1 3 23 24 41 40
+1 3 24 25 42 41
+1 3 25 26 43 42
+1 3 26 27 44 43
+1 3 27 28 45 44
+1 3 28 29 46 45
+1 3 29 30 47 46
+1 3 30 31 48 47
+1 3 31 32 49 48
+1 3 32 33 50 49
+1 3 34 35 52 51
+1 3 35 36 53 52
+1 3 36 37 54 53
+1 3 37 38 55 54
+1 3 38 39 56 55
+1 3 39 40 57 56
+1 3 40 41 58 57
+1 3 41 42 59 58
+1 3 42 43 60 59
+1 3 43 44 61 60
+1 3 44 45 62 61
+1 3 45 46 63 62
+1 3 46 47 64 63
+1 3 47 48 65 64
+1 3 48 49 66 65
+1 3 49 50 67 66
+1 3 51 52 69 68
+1 3 52 53 70 69
+1 3 53 54 71 70
+1 3 54 55 72 71
+1 3 55 56 73 72
+1 3 56 57 74 73
+1 3 57 58 75 74
+1 3 58 59 76 75
+1 3 59 60 77 76
+1 3 60 61 78 77
+1 3 61 62 79 78
+1 3 62 63 80 79
+1 3 63 64 81 80
+1 3 64 65 82 81
+1 3 65 66 83 82
+1 3 66 67 84 83
+
+boundary
+40
+1 1 0 1
+1 1 1 2
+1 1 2 3
+1 1 3 4
+1 1 4 5
+1 1 5 6
+1 1 6 7
+1 1 7 8
+1 1 8 9
+1 1 9 10
+1 1 10 11
+1 1 11 12
+1 1 12 13
+1 1 13 14
+1 1 14 15
+1 1 15 16
+3 1 69 68
+3 1 70 69
+3 1 71 70
+3 1 72 71
+3 1 73 72
+3 1 74 73
+3 1 75 74
+3 1 76 75
+3 1 77 76
+3 1 78 77
+3 1 79 78
+3 1 80 79
+3 1 81 80
+3 1 82 81
+3 1 83 82
+3 1 84 83
+4 1 17 0
+4 1 34 17
+4 1 51 34
+4 1 68 51
+2 1 16 33
+2 1 33 50
+2 1 50 67
+2 1 67 84
+
+vertices
+85
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_2D_P1
+VDim: 2
+Ordering: 1
+
+0 0
+0.25 0
+0 0.25
+0.25 0.25
+0.25 0
+0.5 0
+0.25 0.25
+0.5 0.25
+0.5 0
+0.75 0
+0.5 0.25
+0.75 0.25
+0.75 0
+1 0
+0.75 0.25
+1 0.25
+1 0
+1.25 0
+1 0.25
+1.25 0.25
+1.25 0
+1.5 0
+1.25 0.25
+1.5 0.25
+1.5 0
+1.75 0
+1.5 0.25
+1.75 0.25
+1.75 0
+2 0
+1.75 0.25
+2 0.25
+2 0
+2.25 0
+2 0.25
+2.25 0.25
+2.25 0
+2.5 0
+2.25 0.25
+2.5 0.25
+2.5 0
+2.75 0
+2.5 0.25
+2.75 0.25
+2.75 0
+3 0
+2.75 0.25
+3 0.25
+3 0
+3.25 0
+3 0.25
+3.25 0.25
+3.25 0
+3.5 0
+3.25 0.25
+3.5 0.25
+3.5 0
+3.75 0
+3.5 0.25
+3.75 0.25
+3.75 0
+4 0
+3.75 0.25
+4 0.25
+0 0.25
+0.25 0.25
+0 0.5
+0.25 0.5
+0.25 0.25
+0.5 0.25
+0.25 0.5
+0.5 0.5
+0.5 0.25
+0.75 0.25
+0.5 0.5
+0.75 0.5
+0.75 0.25
+1 0.25
+0.75 0.5
+1 0.5
+1 0.25
+1.25 0.25
+1 0.5
+1.25 0.5
+1.25 0.25
+1.5 0.25
+1.25 0.5
+1.5 0.5
+1.5 0.25
+1.75 0.25
+1.5 0.5
+1.75 0.5
+1.75 0.25
+2 0.25
+1.75 0.5
+2 0.5
+2 0.25
+2.25 0.25
+2 0.5
+2.25 0.5
+2.25 0.25
+2.5 0.25
+2.25 0.5
+2.5 0.5
+2.5 0.25
+2.75 0.25
+2.5 0.5
+2.75 0.5
+2.75 0.25
+3 0.25
+2.75 0.5
+3 0.5
+3 0.25
+3.25 0.25
+3 0.5
+3.25 0.5
+3.25 0.25
+3.5 0.25
+3.25 0.5
+3.5 0.5
+3.5 0.25
+3.75 0.25
+3.5 0.5
+3.75 0.5
+3.75 0.25
+4 0.25
+3.75 0.5
+4 0.5
+0 0.5
+0.25 0.5
+0 0.75
+0.25 0.75
+0.25 0.5
+0.5 0.5
+0.25 0.75
+0.5 0.75
+0.5 0.5
+0.75 0.5
+0.5 0.75
+0.75 0.75
+0.75 0.5
+1 0.5
+0.75 0.75
+1 0.75
+1 0.5
+1.25 0.5
+1 0.75
+1.25 0.75
+1.25 0.5
+1.5 0.5
+1.25 0.75
+1.5 0.75
+1.5 0.5
+1.75 0.5
+1.5 0.75
+1.75 0.75
+1.75 0.5
+2 0.5
+1.75 0.75
+2 0.75
+2 0.5
+2.25 0.5
+2 0.75
+2.25 0.75
+2.25 0.5
+2.5 0.5
+2.25 0.75
+2.5 0.75
+2.5 0.5
+2.75 0.5
+2.5 0.75
+2.75 0.75
+2.75 0.5
+3 0.5
+2.75 0.75
+3 0.75
+3 0.5
+3.25 0.5
+3 0.75
+3.25 0.75
+3.25 0.5
+3.5 0.5
+3.25 0.75
+3.5 0.75
+3.5 0.5
+3.75 0.5
+3.5 0.75
+3.75 0.75
+3.75 0.5
+4 0.5
+3.75 0.75
+4 0.75
+0 0.75
+0.25 0.75
+0 1
+0.25 1
+0.25 0.75
+0.5 0.75
+0.25 1
+0.5 1
+0.5 0.75
+0.75 0.75
+0.5 1
+0.75 1
+0.75 0.75
+1 0.75
+0.75 1
+1 1
+1 0.75
+1.25 0.75
+1 1
+1.25 1
+1.25 0.75
+1.5 0.75
+1.25 1
+1.5 1
+1.5 0.75
+1.75 0.75
+1.5 1
+1.75 1
+1.75 0.75
+2 0.75
+1.75 1
+2 1
+2 0.75
+2.25 0.75
+2 1
+2.25 1
+2.25 0.75
+2.5 0.75
+2.25 1
+2.5 1
+2.5 0.75
+2.75 0.75
+2.5 1
+2.75 1
+2.75 0.75
+3 0.75
+2.75 1
+3 1
+3 0.75
+3.25 0.75
+3 1
+3.25 1
+3.25 0.75
+3.5 0.75
+3.25 1
+3.5 1
+3.5 0.75
+3.75 0.75
+3.5 1
+3.75 1
+3.75 0.75
+4 0.75
+3.75 1
+4 1
diff --git a/examples/linelast/meshes/rod2D_V.mesh b/examples/linelast/meshes/rod2D_V.mesh
new file mode 100644
index 00000000..073f5190
--- /dev/null
+++ b/examples/linelast/meshes/rod2D_V.mesh
@@ -0,0 +1,393 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+2
+
+elements
+64
+1 3 0 1 6 5
+1 3 1 2 7 6
+1 3 2 3 8 7
+1 3 3 4 9 8
+1 3 5 6 11 10
+1 3 6 7 12 11
+1 3 7 8 13 12
+1 3 8 9 14 13
+1 3 10 11 16 15
+1 3 11 12 17 16
+1 3 12 13 18 17
+1 3 13 14 19 18
+1 3 15 16 21 20
+1 3 16 17 22 21
+1 3 17 18 23 22
+1 3 18 19 24 23
+1 3 20 21 26 25
+1 3 21 22 27 26
+1 3 22 23 28 27
+1 3 23 24 29 28
+1 3 25 26 31 30
+1 3 26 27 32 31
+1 3 27 28 33 32
+1 3 28 29 34 33
+1 3 30 31 36 35
+1 3 31 32 37 36
+1 3 32 33 38 37
+1 3 33 34 39 38
+1 3 35 36 41 40
+1 3 36 37 42 41
+1 3 37 38 43 42
+1 3 38 39 44 43
+1 3 40 41 46 45
+1 3 41 42 47 46
+1 3 42 43 48 47
+1 3 43 44 49 48
+1 3 45 46 51 50
+1 3 46 47 52 51
+1 3 47 48 53 52
+1 3 48 49 54 53
+1 3 50 51 56 55
+1 3 51 52 57 56
+1 3 52 53 58 57
+1 3 53 54 59 58
+1 3 55 56 61 60
+1 3 56 57 62 61
+1 3 57 58 63 62
+1 3 58 59 64 63
+1 3 60 61 66 65
+1 3 61 62 67 66
+1 3 62 63 68 67
+1 3 63 64 69 68
+1 3 65 66 71 70
+1 3 66 67 72 71
+1 3 67 68 73 72
+1 3 68 69 74 73
+1 3 70 71 76 75
+1 3 71 72 77 76
+1 3 72 73 78 77
+1 3 73 74 79 78
+1 3 75 76 81 80
+1 3 76 77 82 81
+1 3 77 78 83 82
+1 3 78 79 84 83
+
+boundary
+40
+1 1 0 1
+1 1 1 2
+1 1 2 3
+1 1 3 4
+3 1 81 80
+3 1 82 81
+3 1 83 82
+3 1 84 83
+4 1 5 0
+4 1 10 5
+4 1 15 10
+4 1 20 15
+4 1 25 20
+4 1 30 25
+4 1 35 30
+4 1 40 35
+4 1 45 40
+4 1 50 45
+4 1 55 50
+4 1 60 55
+4 1 65 60
+4 1 70 65
+4 1 75 70
+4 1 80 75
+2 1 4 9
+2 1 9 14
+2 1 14 19
+2 1 19 24
+2 1 24 29
+2 1 29 34
+2 1 34 39
+2 1 39 44
+2 1 44 49
+2 1 49 54
+2 1 54 59
+2 1 59 64
+2 1 64 69
+2 1 69 74
+2 1 74 79
+2 1 79 84
+
+vertices
+85
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_2D_P1
+VDim: 2
+Ordering: 1
+
+0 0
+0.25 0
+0 0.25
+0.25 0.25
+0.25 0
+0.5 0
+0.25 0.25
+0.5 0.25
+0.5 0
+0.75 0
+0.5 0.25
+0.75 0.25
+0.75 0
+1 0
+0.75 0.25
+1 0.25
+0 0.25
+0.25 0.25
+0 0.5
+0.25 0.5
+0.25 0.25
+0.5 0.25
+0.25 0.5
+0.5 0.5
+0.5 0.25
+0.75 0.25
+0.5 0.5
+0.75 0.5
+0.75 0.25
+1 0.25
+0.75 0.5
+1 0.5
+0 0.5
+0.25 0.5
+0 0.75
+0.25 0.75
+0.25 0.5
+0.5 0.5
+0.25 0.75
+0.5 0.75
+0.5 0.5
+0.75 0.5
+0.5 0.75
+0.75 0.75
+0.75 0.5
+1 0.5
+0.75 0.75
+1 0.75
+0 0.75
+0.25 0.75
+0 1
+0.25 1
+0.25 0.75
+0.5 0.75
+0.25 1
+0.5 1
+0.5 0.75
+0.75 0.75
+0.5 1
+0.75 1
+0.75 0.75
+1 0.75
+0.75 1
+1 1
+0 1
+0.25 1
+0 1.25
+0.25 1.25
+0.25 1
+0.5 1
+0.25 1.25
+0.5 1.25
+0.5 1
+0.75 1
+0.5 1.25
+0.75 1.25
+0.75 1
+1 1
+0.75 1.25
+1 1.25
+0 1.25
+0.25 1.25
+0 1.5
+0.25 1.5
+0.25 1.25
+0.5 1.25
+0.25 1.5
+0.5 1.5
+0.5 1.25
+0.75 1.25
+0.5 1.5
+0.75 1.5
+0.75 1.25
+1 1.25
+0.75 1.5
+1 1.5
+0 1.5
+0.25 1.5
+0 1.75
+0.25 1.75
+0.25 1.5
+0.5 1.5
+0.25 1.75
+0.5 1.75
+0.5 1.5
+0.75 1.5
+0.5 1.75
+0.75 1.75
+0.75 1.5
+1 1.5
+0.75 1.75
+1 1.75
+0 1.75
+0.25 1.75
+0 2
+0.25 2
+0.25 1.75
+0.5 1.75
+0.25 2
+0.5 2
+0.5 1.75
+0.75 1.75
+0.5 2
+0.75 2
+0.75 1.75
+1 1.75
+0.75 2
+1 2
+0 2
+0.25 2
+0 2.25
+0.25 2.25
+0.25 2
+0.5 2
+0.25 2.25
+0.5 2.25
+0.5 2
+0.75 2
+0.5 2.25
+0.75 2.25
+0.75 2
+1 2
+0.75 2.25
+1 2.25
+0 2.25
+0.25 2.25
+0 2.5
+0.25 2.5
+0.25 2.25
+0.5 2.25
+0.25 2.5
+0.5 2.5
+0.5 2.25
+0.75 2.25
+0.5 2.5
+0.75 2.5
+0.75 2.25
+1 2.25
+0.75 2.5
+1 2.5
+0 2.5
+0.25 2.5
+0 2.75
+0.25 2.75
+0.25 2.5
+0.5 2.5
+0.25 2.75
+0.5 2.75
+0.5 2.5
+0.75 2.5
+0.5 2.75
+0.75 2.75
+0.75 2.5
+1 2.5
+0.75 2.75
+1 2.75
+0 2.75
+0.25 2.75
+0 3
+0.25 3
+0.25 2.75
+0.5 2.75
+0.25 3
+0.5 3
+0.5 2.75
+0.75 2.75
+0.5 3
+0.75 3
+0.75 2.75
+1 2.75
+0.75 3
+1 3
+0 3
+0.25 3
+0 3.25
+0.25 3.25
+0.25 3
+0.5 3
+0.25 3.25
+0.5 3.25
+0.5 3
+0.75 3
+0.5 3.25
+0.75 3.25
+0.75 3
+1 3
+0.75 3.25
+1 3.25
+0 3.25
+0.25 3.25
+0 3.5
+0.25 3.5
+0.25 3.25
+0.5 3.25
+0.25 3.5
+0.5 3.5
+0.5 3.25
+0.75 3.25
+0.5 3.5
+0.75 3.5
+0.75 3.25
+1 3.25
+0.75 3.5
+1 3.5
+0 3.5
+0.25 3.5
+0 3.75
+0.25 3.75
+0.25 3.5
+0.5 3.5
+0.25 3.75
+0.5 3.75
+0.5 3.5
+0.75 3.5
+0.5 3.75
+0.75 3.75
+0.75 3.5
+1 3.5
+0.75 3.75
+1 3.75
+0 3.75
+0.25 3.75
+0 4
+0.25 4
+0.25 3.75
+0.5 3.75
+0.25 4
+0.5 4
+0.5 3.75
+0.75 3.75
+0.5 4
+0.75 4
+0.75 3.75
+1 3.75
+0.75 4
+1 4
diff --git a/include/linelast_solver.hpp b/include/linelast_solver.hpp
index 88d089d5..0ae4b356 100644
--- a/include/linelast_solver.hpp
+++ b/include/linelast_solver.hpp
@@ -50,7 +50,7 @@ class LinElastSolver : public MultiBlockSolver
    double kappa = -1.0;
 
    // Initial positions
-   VectorFunctionCoefficient *init_x;
+   VectorFunctionCoefficient *init_x = NULL;
 
 public:
    LinElastSolver();
diff --git a/include/parameterized_problem.hpp b/include/parameterized_problem.hpp
index b25aefd6..c9d80238 100644
--- a/include/parameterized_problem.hpp
+++ b/include/parameterized_problem.hpp
@@ -228,6 +228,18 @@ class LinElastDisp : public LinElastProblem
    LinElastDisp();
 };
 
+class LinElastDispLCantilever : public LinElastProblem
+{
+public:
+   LinElastDispLCantilever();
+};
+
+class LinElastDispLattice : public LinElastProblem
+{
+public:
+   LinElastDispLattice();
+};
+
 ParameterizedProblem* InitParameterizedProblem();
 
 #endif
diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp
index 2512900c..d90a6fef 100644
--- a/src/linelast_solver.cpp
+++ b/src/linelast_solver.cpp
@@ -56,6 +56,7 @@ LinElastSolver::~LinElastSolver()
    delete globalMat_hypre;
    delete mumps;
    delete init_x;
+   
 }
 
 void LinElastSolver::SetupIC(std::function<void(const Vector &, Vector &)> F)
@@ -524,8 +525,8 @@ void LinElastSolver::BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp)
          bdr_marker = 0;
          bdr_marker[comp->bdr_attributes[b] - 1] = 1;
          BilinearForm a_comp(fes_comp[c]);
-         a_comp.AddBdrFaceIntegrator(new DGElasticityIntegrator(*(lambda_c[c]), *(mu_c[c]), alpha, kappa), *(bdr_markers[b]));
 
+         a_comp.AddBdrFaceIntegrator(new DGElasticityIntegrator(*(lambda_c[c]), *(mu_c[c]), alpha, kappa), bdr_marker);
          a_comp.Assemble();
          a_comp.Finalize();
 
diff --git a/src/parameterized_problem.cpp b/src/parameterized_problem.cpp
index cd9e4cc8..9d0ca6f6 100644
--- a/src/parameterized_problem.cpp
+++ b/src/parameterized_problem.cpp
@@ -179,6 +179,11 @@ void init_disp(const Vector &x, Vector &u)
    u(u.Size()-1) = -0.2*x(0)*rdisp_f;
 }
 
+void init_disp_lcantilever(const Vector &x, Vector &u)
+{
+   u = 0.0;
+   u(u.Size()-1) = -0.2*(x(u.Size()-1) - 5.0)*rdisp_f;
+}
 
 }  // namespace linelast_disp
 
@@ -274,6 +279,14 @@ ParameterizedProblem* InitParameterizedProblem()
    {
       problem = new LinElastDisp();
    }
+   else if (problem_name == "linelast_disp_lcantilever")
+   {
+      problem = new LinElastDispLCantilever();
+   }
+   else if (problem_name == "linelast_disp_lattice")
+   {
+      problem = new LinElastDispLattice();
+   }
    else
    {
       mfem_error("Unknown parameterized problem name!\n");
@@ -597,3 +610,78 @@ LinElastDisp::LinElastDisp()
    param_ptr[1] = &(function_factory::linelast_problem::_lambda);
    param_ptr[2] = &(function_factory::linelast_problem::_mu);
 }
+
+LinElastDispLCantilever::LinElastDispLCantilever()
+    : LinElastProblem()
+{
+   // pointer to static function.
+   bdr_type.SetSize(2);
+   battr.SetSize(2);
+   vector_bdr_ptr.SetSize(2);
+   for (size_t i = 0; i < vector_bdr_ptr.Size(); i++)
+   {
+   bdr_type[i] = LinElastProblem::DIRICHLET;
+   battr[i] = i+1;
+   vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp_lcantilever);
+   }
+   
+   // Set materials
+   general_scalar_ptr.SetSize(2);
+   general_scalar_ptr[0] = function_factory::linelast_problem::lambda;
+   general_scalar_ptr[1] = function_factory::linelast_problem::mu;
+   
+   // Default values.
+   function_factory::linelast_disp::rdisp_f = 1.0;
+   function_factory::linelast_problem::_lambda = 1.0;
+   function_factory::linelast_problem::_mu = 1.0;
+
+   param_map["rdisp_f"] = 0;
+   param_map["lambda"] = 1;
+   param_map["mu"] = 2;
+
+   param_ptr.SetSize(3);
+   param_ptr[0] = &(function_factory::linelast_disp::rdisp_f);
+   param_ptr[1] = &(function_factory::linelast_problem::_lambda);
+   param_ptr[2] = &(function_factory::linelast_problem::_mu);
+
+   general_vector_ptr.SetSize(1);
+   general_vector_ptr[0] = NULL;
+}
+
+LinElastDispLattice::LinElastDispLattice()
+    : LinElastProblem()
+{
+   // pointer to static function.
+   bdr_type.SetSize(2);
+   battr.SetSize(2);
+   vector_bdr_ptr.SetSize(2);
+   for (size_t i = 0; i < vector_bdr_ptr.Size(); i++)
+   {
+   bdr_type[i] = LinElastProblem::DIRICHLET;
+   battr[i] = i+1;
+   vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp);
+   }
+   
+   // Set materials
+   general_scalar_ptr.SetSize(2);
+   general_scalar_ptr[0] = function_factory::linelast_problem::lambda;
+   general_scalar_ptr[1] = function_factory::linelast_problem::mu;
+
+   // Default values.
+   function_factory::linelast_disp::rdisp_f = 1.0;
+   function_factory::linelast_problem::_lambda = 1.0;
+   function_factory::linelast_problem::_mu = 1.0;
+
+   param_map["rdisp_f"] = 0;
+   param_map["lambda"] = 1;
+   param_map["mu"] = 2;
+
+   param_ptr.SetSize(3);
+   param_ptr[0] = &(function_factory::linelast_disp::rdisp_f);
+   param_ptr[1] = &(function_factory::linelast_problem::_lambda);
+   param_ptr[2] = &(function_factory::linelast_problem::_mu);
+
+   general_vector_ptr.SetSize(1);
+   general_vector_ptr[0] = NULL;
+
+}
diff --git a/utils/python/linelast_comp.py b/utils/python/linelast_comp.py
new file mode 100644
index 00000000..0a54b4f5
--- /dev/null
+++ b/utils/python/linelast_comp.py
@@ -0,0 +1,226 @@
+# Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details.
+#
+# SPDX-License-Identifier: MIT
+
+import numpy as np
+import h5py
+
+def SimpleL():
+    n_mesh = 3
+    mesh_type = [0, 1, 2]
+    mesh_configs = np.zeros([n_mesh, 6])
+    mesh_configs[1,:] = [1., 0., 0., 0., 0., 0.]
+    mesh_configs[2,:] = [0., 1., 0., 0., 0., 0.]
+
+    # interface data
+    # mesh1 / mesh2 / battr1 / battr2 / port_idx
+    if_data = np.zeros([2, 5])
+    if_data[0, :] = [0, 1, 2, 4, 0]
+    if_data[1, :] = [0, 2, 3, 1, 1]
+
+    # boundary attributes
+    # global_battr / mesh_idx / comp_battr
+    bdr_data = []
+
+    # applied loads and displacements
+    bdr_data += [[1, 1, 2]]
+    bdr_data += [[2, 2, 3]]
+
+    # homogenous neumann
+    bdr_data += [[3, 0, 4]]
+    bdr_data += [[3, 0, 1]]
+    bdr_data += [[3, 1, 1]]
+    bdr_data += [[3, 1, 3]]
+    bdr_data += [[3, 2, 2]]
+    bdr_data += [[3, 2, 4]]
+
+    bdr_data = np.array(bdr_data)
+    print(bdr_data.shape)
+
+    filename = "linelast.simpleL.h5"
+    with h5py.File(filename, 'w') as f:
+        # c++ currently cannot read datasets of string.
+        # change to multiple attributes, only as a temporary implementation.
+        grp = f.create_group("components")
+        grp.attrs["number_of_components"] = 3
+        grp.attrs["0"] = "joint2D"
+        grp.attrs["1"] = "rod2D_H"
+        grp.attrs["2"] = "rod2D_V"
+        # component index of each mesh
+        grp.create_dataset("meshes", (n_mesh,), data=mesh_type)
+        # 3-dimension vector for translation / rotation
+        grp.create_dataset("configuration", mesh_configs.shape, data=mesh_configs)
+
+        grp = f.create_group("ports")
+        grp.attrs["number_of_references"] = 2
+        grp.attrs["0"] = "port1"
+        grp.attrs["1"] = "port2"
+        grp.create_dataset("interface", if_data.shape, data=if_data)
+
+        port = grp.create_group("port1")
+        port.attrs["comp1"] = "joint2D"
+        port.attrs["comp2"] = "rod2D_H"
+        port.attrs["attr1"] = 2
+        port.attrs["attr2"] = 4
+        port.create_dataset("comp2_configuration", (6,), data=[1., 0., 0., 0., 0., 0.])
+
+        port = grp.create_group("port2")
+        port.attrs["comp1"] = "joint2D"
+        port.attrs["comp2"] = "rod2D_V"
+        port.attrs["attr1"] = 3
+        port.attrs["attr2"] = 1
+        port.create_dataset("comp2_configuration", (6,), data=[0., 1., 0., 0., 0., 0.])
+
+        # boundary attributes
+        f.create_dataset("boundary", bdr_data.shape, data=bdr_data)
+    return
+
+def grid_mesh_configs(nx, ny, l, w):
+    bdr_data = []
+    if_data = []
+    n_joint = (nx + 1) * (ny + 1)
+    n_rod_H = (ny + 1) * nx
+    n_rod_V = (nx + 1) * ny
+    n_rod = n_rod_H + n_rod_V
+    n_mesh = n_joint + n_rod
+    mesh_configs = np.zeros([n_mesh, 6])
+    mesh_type = np.concatenate((np.full((int(n_joint),), 0.0),np.full((int(n_rod_H),), 1.0),np.full((int(n_rod_V),), 2.0)))
+    # Setup joints
+    for i in range(n_joint):
+        xi = (i + nx+1) % (nx+1)
+        yi = np.floor(i / (nx+1))
+        mesh_configs[i,:] = [xi*(l+w), yi*(l+w), 0., 0., 0., 0.]
+
+        # Boundary check
+        # global_battr / mesh_idx / comp_battr
+        if xi == 0.0:
+            bdr_data += [[1, i, 4]] # constrain one end
+        elif xi == nx:
+            bdr_data += [[2, i, 2]] # constrain one end
+        if yi==0.0:
+            bdr_data += [[3, i, 1]] # free boundary
+        elif yi==ny:
+            bdr_data += [[3, i, 3]] # free boundary
+        
+        # Interface check
+        # mesh1 / mesh2 / battr1 / battr2 / port_idx
+        # Case port idx = 0
+        if xi < nx:
+            if_data += [[i, i + n_joint - yi, 2, 4, 0]]
+        # Case port idx = 1
+        if xi > 0:
+            if_data += [[i, i + n_joint - yi - 1, 4, 2, 1]]
+        # Case port idx = 2
+        if yi < ny:
+            if_data += [[i, n_joint + n_rod_H + i, 3, 1, 2]]
+        # Case port idx = 3
+        if yi > 0:
+            if_data += [[i, n_joint + n_rod_H + i - (nx+1), 1, 3, 3]]
+
+    # Setup horizontal rods
+    for i in range(n_rod_H):
+        xi = (i + nx) % (nx)
+        yi = np.floor(i / (nx))
+        mesh_configs[n_joint + i,:] = [xi*(l+w) + w, yi*(l+w), 0., 0., 0., 0.]
+        # Boundary check
+        bdr_data += [[3, n_joint + i, 1]] # free boundary
+        bdr_data += [[3, n_joint + i, 3]] # free boundary
+    
+    # Setup vertical rods
+    for i in range(n_rod_V):
+        xi = (i + nx+1) % (nx+1)
+        yi = np.floor(i / (nx+1))
+        mesh_configs[n_joint + n_rod_H + i,:] = [xi*(l+w), yi*(l+w) + w, 0., 0., 0., 0.]
+
+        # Boundary check
+        if xi == 0.0:
+            bdr_data += [[1, n_joint + n_rod_H + i, 4]] # constrain one end
+            bdr_data += [[3, n_joint + n_rod_H + i, 2]] # free boundary
+        elif xi == nx:
+            bdr_data += [[2, n_joint + n_rod_H + i, 2]] # constrain one end
+            bdr_data += [[3, n_joint + n_rod_H + i, 4]] # free boundary
+        else:
+            bdr_data += [[3, n_joint + n_rod_H + i, 2]] # constrain one end
+            bdr_data += [[3, n_joint + n_rod_H + i, 4]] # free boundary
+
+    return mesh_configs, bdr_data, if_data, mesh_type, n_mesh
+
+
+def LatticeCantilever(nx, ny):
+    # nx and ny are the number of sections in x- and y-direction, respectively
+    l = 4.0
+    w = 1.0
+    mesh_configs, bdr_data, if_data, mesh_type, n_mesh = grid_mesh_configs(nx, ny, l, w)
+
+    # interface data
+    if_data = np.array(if_data)
+
+    # boundary attributes
+    bdr_data = np.array(bdr_data)
+
+    filename = "linelast.lattice.h5"
+    with h5py.File(filename, 'w') as f:
+        # c++ currently cannot read datasets of string.
+        # change to multiple attributes, only as a temporary implementation.
+        grp = f.create_group("components")
+        grp.attrs["number_of_components"] = 3
+        grp.attrs["0"] = "joint2D"
+        grp.attrs["1"] = "rod2D_H"
+        grp.attrs["2"] = "rod2D_V"
+        # component index of each mesh
+        grp.create_dataset("meshes", (n_mesh,), data=mesh_type)
+        # 3-dimension vector for translation / rotation
+        grp.create_dataset("configuration", mesh_configs.shape, data=mesh_configs)
+
+        grp = f.create_group("ports")
+        grp.attrs["number_of_references"] = 4
+        grp.attrs["0"] = "port1"
+        grp.attrs["1"] = "port2"
+        grp.attrs["2"] = "port3"
+        grp.attrs["3"] = "port4"
+        grp.create_dataset("interface", if_data.shape, data=if_data)
+
+        port = grp.create_group("port1")
+        port.attrs["comp1"] = "joint2D"
+        port.attrs["comp2"] = "rod2D_H"
+        port.attrs["attr1"] = 2
+        port.attrs["attr2"] = 4
+        port.create_dataset("comp2_configuration", (6,), data=[w, 0., 0., 0., 0., 0.])
+
+        port = grp.create_group("port2")
+        port.attrs["comp1"] = "joint2D"
+        port.attrs["comp2"] = "rod2D_H"
+        port.attrs["attr1"] = 4
+        port.attrs["attr2"] = 2
+        port.create_dataset("comp2_configuration", (6,), data=[-l, 0., 0., 0., 0., 0.])
+
+        port = grp.create_group("port3")
+        port.attrs["comp1"] = "joint2D"
+        port.attrs["comp2"] = "rod2D_V"
+        port.attrs["attr1"] = 3
+        port.attrs["attr2"] = 1
+        port.create_dataset("comp2_configuration", (6,), data=[0., w, 0., 0., 0., 0.])
+
+        port = grp.create_group("port4")
+        port.attrs["comp1"] = "joint2D"
+        port.attrs["comp2"] = "rod2D_V"
+        port.attrs["attr1"] = 1
+        port.attrs["attr2"] = 3
+        port.create_dataset("comp2_configuration", (6,), data=[0., -l, 0., 0., 0., 0.])
+
+        # boundary attributes
+        f.create_dataset("boundary", bdr_data.shape, data=bdr_data)
+    return
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) == 1:
+        SimpleL()
+    else:
+        name = sys.argv[1]
+        if name == "simple_l":
+            SimpleL()
+        elif name == "lattice_cantilever":
+            nx = int(sys.argv[2])
+            ny = int(sys.argv[3])
+            LatticeCantilever(nx, ny)