From e3cf21faabc34ec159aaa7eb7490c00b7b17767a Mon Sep 17 00:00:00 2001 From: fatima Date: Tue, 4 Aug 2026 16:02:57 +0200 Subject: [PATCH] Add 3D beam 3-point flexion validation: SOFA vs FreeFEM --- .../validation/3D/3D_pointFlex/beam3d_tet.msh | 534 ++++++++++++++++++ .../3D_pointFlex/comparaison_script3d_3pt.py | 189 +++++++ .../3D_pointFlex/freefem_beam3d_threept.edp | 120 ++++ .../3D_pointFlex/params_beam3d_threept.json | 5 + .../3D/3D_pointFlex/sofa_beam3d_threept.py | 192 +++++++ 5 files changed, 1040 insertions(+) create mode 100644 examples/Freefem/validation/3D/3D_pointFlex/beam3d_tet.msh create mode 100644 examples/Freefem/validation/3D/3D_pointFlex/comparaison_script3d_3pt.py create mode 100644 examples/Freefem/validation/3D/3D_pointFlex/freefem_beam3d_threept.edp create mode 100644 examples/Freefem/validation/3D/3D_pointFlex/params_beam3d_threept.json create mode 100644 examples/Freefem/validation/3D/3D_pointFlex/sofa_beam3d_threept.py diff --git a/examples/Freefem/validation/3D/3D_pointFlex/beam3d_tet.msh b/examples/Freefem/validation/3D/3D_pointFlex/beam3d_tet.msh new file mode 100644 index 00000000..264955ad --- /dev/null +++ b/examples/Freefem/validation/3D/3D_pointFlex/beam3d_tet.msh @@ -0,0 +1,534 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +7 +2 1 "Fixed" +2 2 "Bottom" +2 3 "Right" +2 4 "Top" +2 5 "Front" +2 6 "Back" +3 7 "Beam" +$EndPhysicalNames +$Nodes +99 +1 0 0 0 +2 1 0 0 +3 1 0.2 0 +4 0 0.2 0 +5 0 0 0.2 +6 1 0 0.2 +7 1 0.2 0.2 +8 0 0.2 0.2 +9 0.09999999999981413 0 0 +10 0.1999999999995569 0 0 +11 0.299999999999265 0 0 +12 0.3999999999989731 0 0 +13 0.4999999999986921 0 0 +14 0.599999999998945 0 0 +15 0.6999999999992088 0 0 +16 0.7999999999994725 0 0 +17 0.8999999999997362 0 0 +18 1 0.0999999999997371 0 +19 0.8999999999995836 0.2 0 +20 0.8 0.2 0 +21 0.7000000000006938 0.2 0 +22 0.6000000000013874 0.2 0 +23 0.5000000000020595 0.2 0 +24 0.400000000001665 0.2 0 +25 0.3000000000012487 0.2 0 +26 0.2000000000008325 0.2 0 +27 0.1000000000004162 0.2 0 +28 0 0.100000000000274 0 +29 0.09999999999981413 0 0.2 +30 0.1999999999995569 0 0.2 +31 0.299999999999265 0 0.2 +32 0.3999999999989731 0 0.2 +33 0.4999999999986921 0 0.2 +34 0.599999999998945 0 0.2 +35 0.6999999999992088 0 0.2 +36 0.7999999999994725 0 0.2 +37 0.8999999999997362 0 0.2 +38 1 0.0999999999997371 0.2 +39 0.8999999999995836 0.2 0.2 +40 0.8 0.2 0.2 +41 0.7000000000006938 0.2 0.2 +42 0.6000000000013874 0.2 0.2 +43 0.5000000000020595 0.2 0.2 +44 0.400000000001665 0.2 0.2 +45 0.3000000000012487 0.2 0.2 +46 0.2000000000008325 0.2 0.2 +47 0.1000000000004162 0.2 0.2 +48 0 0.100000000000274 0.2 +49 0 0 0.0999999999997371 +50 1 0 0.0999999999997371 +51 1 0.2 0.0999999999997371 +52 0 0.2 0.0999999999997371 +53 0.1000000000001152 0.1000000000002203 0 +54 0.2000000000001947 0.1000000000001666 0 +55 0.3000000000002568 0.1000000000001129 0 +56 0.4000000000003192 0.1000000000000592 0 +57 0.5000000000003757 0.1000000000000055 0 +58 0.6000000000001662 0.09999999999995185 0 +59 0.6999999999999513 0.09999999999989817 0 +60 0.7999999999997363 0.09999999999984449 0 +61 0.89999999999966 0.09999999999979081 0 +62 0.1000000000001152 0.1000000000002203 0.2 +63 0.2000000000001947 0.1000000000001666 0.2 +64 0.3000000000002568 0.1000000000001129 0.2 +65 0.4000000000003192 0.1000000000000592 0.2 +66 0.5000000000003757 0.1000000000000055 0.2 +67 0.6000000000001662 0.09999999999995185 0.2 +68 0.6999999999999513 0.09999999999989817 0.2 +69 0.7999999999997363 0.09999999999984449 0.2 +70 0.89999999999966 0.09999999999979081 0.2 +71 0.09999999999981413 0 0.0999999999997371 +72 0.1999999999995569 0 0.0999999999997371 +73 0.299999999999265 0 0.0999999999997371 +74 0.3999999999989731 0 0.0999999999997371 +75 0.4999999999986922 0 0.0999999999997371 +76 0.599999999998945 0 0.0999999999997371 +77 0.6999999999992088 0 0.0999999999997371 +78 0.7999999999994725 0 0.0999999999997371 +79 0.8999999999997362 0 0.0999999999997371 +80 0.8999999999995836 0.2 0.0999999999997371 +81 0.8 0.2 0.0999999999997371 +82 0.7000000000006938 0.2 0.0999999999997371 +83 0.6000000000013874 0.2 0.0999999999997371 +84 0.5000000000020595 0.2 0.0999999999997371 +85 0.400000000001665 0.2 0.0999999999997371 +86 0.3000000000012487 0.2 0.0999999999997371 +87 0.2000000000008325 0.2 0.0999999999997371 +88 0.1000000000004162 0.2 0.0999999999997371 +89 1 0.0999999999997371 0.0999999999997371 +90 0 0.100000000000274 0.0999999999997371 +91 0.1000000000001152 0.1000000000002202 0.09999999999973708 +92 0.2000000000001948 0.1000000000001666 0.09999999999973708 +93 0.3000000000002568 0.1000000000001129 0.09999999999973709 +94 0.4000000000003192 0.1000000000000592 0.09999999999973706 +95 0.5000000000003757 0.1000000000000055 0.0999999999997371 +96 0.6000000000001662 0.09999999999995182 0.09999999999973708 +97 0.6999999999999513 0.09999999999989814 0.09999999999973708 +98 0.799999999999736 0.09999999999984448 0.09999999999973708 +99 0.8999999999996604 0.09999999999979081 0.09999999999973708 +$EndNodes +$Elements +416 +1 2 2 5 1 1 9 28 +2 2 2 5 1 28 9 53 +3 2 2 5 1 28 53 4 +4 2 2 5 1 4 53 27 +5 2 2 5 1 9 10 53 +6 2 2 5 1 53 10 54 +7 2 2 5 1 53 54 27 +8 2 2 5 1 27 54 26 +9 2 2 5 1 10 11 54 +10 2 2 5 1 54 11 55 +11 2 2 5 1 54 55 26 +12 2 2 5 1 26 55 25 +13 2 2 5 1 11 12 55 +14 2 2 5 1 55 12 56 +15 2 2 5 1 55 56 25 +16 2 2 5 1 25 56 24 +17 2 2 5 1 12 13 56 +18 2 2 5 1 56 13 57 +19 2 2 5 1 56 57 24 +20 2 2 5 1 24 57 23 +21 2 2 5 1 13 14 57 +22 2 2 5 1 57 14 58 +23 2 2 5 1 57 58 23 +24 2 2 5 1 23 58 22 +25 2 2 5 1 14 15 58 +26 2 2 5 1 58 15 59 +27 2 2 5 1 58 59 22 +28 2 2 5 1 22 59 21 +29 2 2 5 1 15 16 59 +30 2 2 5 1 59 16 60 +31 2 2 5 1 59 60 21 +32 2 2 5 1 21 60 20 +33 2 2 5 1 16 17 60 +34 2 2 5 1 60 17 61 +35 2 2 5 1 60 61 20 +36 2 2 5 1 20 61 19 +37 2 2 5 1 17 2 61 +38 2 2 5 1 61 2 18 +39 2 2 5 1 61 18 19 +40 2 2 5 1 19 18 3 +41 2 2 6 2 5 29 48 +42 2 2 6 2 48 29 62 +43 2 2 6 2 48 62 8 +44 2 2 6 2 8 62 47 +45 2 2 6 2 29 30 62 +46 2 2 6 2 62 30 63 +47 2 2 6 2 62 63 47 +48 2 2 6 2 47 63 46 +49 2 2 6 2 30 31 63 +50 2 2 6 2 63 31 64 +51 2 2 6 2 63 64 46 +52 2 2 6 2 46 64 45 +53 2 2 6 2 31 32 64 +54 2 2 6 2 64 32 65 +55 2 2 6 2 64 65 45 +56 2 2 6 2 45 65 44 +57 2 2 6 2 32 33 65 +58 2 2 6 2 65 33 66 +59 2 2 6 2 65 66 44 +60 2 2 6 2 44 66 43 +61 2 2 6 2 33 34 66 +62 2 2 6 2 66 34 67 +63 2 2 6 2 66 67 43 +64 2 2 6 2 43 67 42 +65 2 2 6 2 34 35 67 +66 2 2 6 2 67 35 68 +67 2 2 6 2 67 68 42 +68 2 2 6 2 42 68 41 +69 2 2 6 2 35 36 68 +70 2 2 6 2 68 36 69 +71 2 2 6 2 68 69 41 +72 2 2 6 2 41 69 40 +73 2 2 6 2 36 37 69 +74 2 2 6 2 69 37 70 +75 2 2 6 2 69 70 40 +76 2 2 6 2 40 70 39 +77 2 2 6 2 37 6 70 +78 2 2 6 2 70 6 38 +79 2 2 6 2 70 38 39 +80 2 2 6 2 39 38 7 +81 2 2 2 3 1 9 49 +82 2 2 2 3 49 9 71 +83 2 2 2 3 49 71 5 +84 2 2 2 3 5 71 29 +85 2 2 2 3 9 10 71 +86 2 2 2 3 71 10 72 +87 2 2 2 3 71 72 29 +88 2 2 2 3 29 72 30 +89 2 2 2 3 10 11 72 +90 2 2 2 3 72 11 73 +91 2 2 2 3 72 73 30 +92 2 2 2 3 30 73 31 +93 2 2 2 3 11 12 73 +94 2 2 2 3 73 12 74 +95 2 2 2 3 73 74 31 +96 2 2 2 3 31 74 32 +97 2 2 2 3 12 13 74 +98 2 2 2 3 74 13 75 +99 2 2 2 3 74 75 32 +100 2 2 2 3 32 75 33 +101 2 2 2 3 13 14 75 +102 2 2 2 3 75 14 76 +103 2 2 2 3 75 76 33 +104 2 2 2 3 33 76 34 +105 2 2 2 3 14 15 76 +106 2 2 2 3 76 15 77 +107 2 2 2 3 76 77 34 +108 2 2 2 3 34 77 35 +109 2 2 2 3 15 16 77 +110 2 2 2 3 77 16 78 +111 2 2 2 3 77 78 35 +112 2 2 2 3 35 78 36 +113 2 2 2 3 16 17 78 +114 2 2 2 3 78 17 79 +115 2 2 2 3 78 79 36 +116 2 2 2 3 36 79 37 +117 2 2 2 3 17 2 79 +118 2 2 2 3 79 2 50 +119 2 2 2 3 79 50 37 +120 2 2 2 3 37 50 6 +121 2 2 3 5 2 18 50 +122 2 2 3 5 50 18 89 +123 2 2 3 5 50 89 6 +124 2 2 3 5 6 89 38 +125 2 2 3 5 18 3 89 +126 2 2 3 5 89 3 51 +127 2 2 3 5 89 51 38 +128 2 2 3 5 38 51 7 +129 4 2 7 1 1 9 28 49 +130 4 2 7 1 9 28 49 71 +131 4 2 7 1 49 71 28 90 +132 4 2 7 1 9 28 71 53 +133 4 2 7 1 28 90 71 53 +134 4 2 7 1 71 90 91 53 +135 4 2 7 1 49 71 90 5 +136 4 2 7 1 71 90 5 29 +137 4 2 7 1 5 29 90 48 +138 4 2 7 1 71 90 29 91 +139 4 2 7 1 90 48 29 91 +140 4 2 7 1 29 48 62 91 +141 4 2 7 1 28 53 4 90 +142 4 2 7 1 53 4 90 91 +143 4 2 7 1 90 91 4 52 +144 4 2 7 1 53 4 91 27 +145 4 2 7 1 4 52 91 27 +146 4 2 7 1 91 52 88 27 +147 4 2 7 1 90 91 52 48 +148 4 2 7 1 91 52 48 62 +149 4 2 7 1 48 62 52 8 +150 4 2 7 1 91 52 62 88 +151 4 2 7 1 52 8 62 88 +152 4 2 7 1 62 8 47 88 +153 4 2 7 1 9 10 53 71 +154 4 2 7 1 10 53 71 72 +155 4 2 7 1 71 72 53 91 +156 4 2 7 1 10 53 72 54 +157 4 2 7 1 53 91 72 54 +158 4 2 7 1 72 91 92 54 +159 4 2 7 1 71 72 91 29 +160 4 2 7 1 72 91 29 30 +161 4 2 7 1 29 30 91 62 +162 4 2 7 1 72 91 30 92 +163 4 2 7 1 91 62 30 92 +164 4 2 7 1 30 62 63 92 +165 4 2 7 1 53 54 27 91 +166 4 2 7 1 54 27 91 92 +167 4 2 7 1 91 92 27 88 +168 4 2 7 1 54 27 92 26 +169 4 2 7 1 27 88 92 26 +170 4 2 7 1 92 88 87 26 +171 4 2 7 1 91 92 88 62 +172 4 2 7 1 92 88 62 63 +173 4 2 7 1 62 63 88 47 +174 4 2 7 1 92 88 63 87 +175 4 2 7 1 88 47 63 87 +176 4 2 7 1 63 47 46 87 +177 4 2 7 1 10 11 54 72 +178 4 2 7 1 11 54 72 73 +179 4 2 7 1 72 73 54 92 +180 4 2 7 1 11 54 73 55 +181 4 2 7 1 54 92 73 55 +182 4 2 7 1 73 92 93 55 +183 4 2 7 1 72 73 92 30 +184 4 2 7 1 73 92 30 31 +185 4 2 7 1 30 31 92 63 +186 4 2 7 1 73 92 31 93 +187 4 2 7 1 92 63 31 93 +188 4 2 7 1 31 63 64 93 +189 4 2 7 1 54 55 26 92 +190 4 2 7 1 55 26 92 93 +191 4 2 7 1 92 93 26 87 +192 4 2 7 1 55 26 93 25 +193 4 2 7 1 26 87 93 25 +194 4 2 7 1 93 87 86 25 +195 4 2 7 1 92 93 87 63 +196 4 2 7 1 93 87 63 64 +197 4 2 7 1 63 64 87 46 +198 4 2 7 1 93 87 64 86 +199 4 2 7 1 87 46 64 86 +200 4 2 7 1 64 46 45 86 +201 4 2 7 1 11 12 55 73 +202 4 2 7 1 12 55 73 74 +203 4 2 7 1 73 74 55 93 +204 4 2 7 1 12 55 74 56 +205 4 2 7 1 55 93 74 56 +206 4 2 7 1 74 93 94 56 +207 4 2 7 1 73 74 93 31 +208 4 2 7 1 74 93 31 32 +209 4 2 7 1 31 32 93 64 +210 4 2 7 1 74 93 32 94 +211 4 2 7 1 93 64 32 94 +212 4 2 7 1 32 64 65 94 +213 4 2 7 1 55 56 25 93 +214 4 2 7 1 56 25 93 94 +215 4 2 7 1 93 94 25 86 +216 4 2 7 1 56 25 94 24 +217 4 2 7 1 25 86 94 24 +218 4 2 7 1 94 86 85 24 +219 4 2 7 1 93 94 86 64 +220 4 2 7 1 94 86 64 65 +221 4 2 7 1 64 65 86 45 +222 4 2 7 1 94 86 65 85 +223 4 2 7 1 86 45 65 85 +224 4 2 7 1 65 45 44 85 +225 4 2 7 1 12 13 56 74 +226 4 2 7 1 13 56 74 75 +227 4 2 7 1 74 75 56 94 +228 4 2 7 1 13 56 75 57 +229 4 2 7 1 56 94 75 57 +230 4 2 7 1 75 94 95 57 +231 4 2 7 1 74 75 94 32 +232 4 2 7 1 75 94 32 33 +233 4 2 7 1 32 33 94 65 +234 4 2 7 1 75 94 33 95 +235 4 2 7 1 94 65 33 95 +236 4 2 7 1 33 65 66 95 +237 4 2 7 1 56 57 24 94 +238 4 2 7 1 57 24 94 95 +239 4 2 7 1 94 95 24 85 +240 4 2 7 1 57 24 95 23 +241 4 2 7 1 24 85 95 23 +242 4 2 7 1 95 85 84 23 +243 4 2 7 1 94 95 85 65 +244 4 2 7 1 95 85 65 66 +245 4 2 7 1 65 66 85 44 +246 4 2 7 1 95 85 66 84 +247 4 2 7 1 85 44 66 84 +248 4 2 7 1 66 44 43 84 +249 4 2 7 1 13 14 57 75 +250 4 2 7 1 14 57 75 76 +251 4 2 7 1 75 76 57 95 +252 4 2 7 1 14 57 76 58 +253 4 2 7 1 57 95 76 58 +254 4 2 7 1 76 95 96 58 +255 4 2 7 1 75 76 95 33 +256 4 2 7 1 76 95 33 34 +257 4 2 7 1 33 34 95 66 +258 4 2 7 1 76 95 34 96 +259 4 2 7 1 95 66 34 96 +260 4 2 7 1 34 66 67 96 +261 4 2 7 1 57 58 23 95 +262 4 2 7 1 58 23 95 96 +263 4 2 7 1 95 96 23 84 +264 4 2 7 1 58 23 96 22 +265 4 2 7 1 23 84 96 22 +266 4 2 7 1 96 84 83 22 +267 4 2 7 1 95 96 84 66 +268 4 2 7 1 96 84 66 67 +269 4 2 7 1 66 67 84 43 +270 4 2 7 1 96 84 67 83 +271 4 2 7 1 84 43 67 83 +272 4 2 7 1 67 43 42 83 +273 4 2 7 1 14 15 58 76 +274 4 2 7 1 15 58 76 77 +275 4 2 7 1 76 77 58 96 +276 4 2 7 1 15 58 77 59 +277 4 2 7 1 58 96 77 59 +278 4 2 7 1 77 96 97 59 +279 4 2 7 1 76 77 96 34 +280 4 2 7 1 77 96 34 35 +281 4 2 7 1 34 35 96 67 +282 4 2 7 1 77 96 35 97 +283 4 2 7 1 96 67 35 97 +284 4 2 7 1 35 67 68 97 +285 4 2 7 1 58 59 22 96 +286 4 2 7 1 59 22 96 97 +287 4 2 7 1 96 97 22 83 +288 4 2 7 1 59 22 97 21 +289 4 2 7 1 22 83 97 21 +290 4 2 7 1 97 83 82 21 +291 4 2 7 1 96 97 83 67 +292 4 2 7 1 97 83 67 68 +293 4 2 7 1 67 68 83 42 +294 4 2 7 1 97 83 68 82 +295 4 2 7 1 83 42 68 82 +296 4 2 7 1 68 42 41 82 +297 4 2 7 1 15 16 59 77 +298 4 2 7 1 16 59 77 78 +299 4 2 7 1 77 78 59 97 +300 4 2 7 1 16 59 78 60 +301 4 2 7 1 59 97 78 60 +302 4 2 7 1 78 97 98 60 +303 4 2 7 1 77 78 97 35 +304 4 2 7 1 78 97 35 36 +305 4 2 7 1 35 36 97 68 +306 4 2 7 1 78 97 36 98 +307 4 2 7 1 97 68 36 98 +308 4 2 7 1 36 68 69 98 +309 4 2 7 1 59 60 21 97 +310 4 2 7 1 60 21 97 98 +311 4 2 7 1 97 98 21 82 +312 4 2 7 1 60 21 98 20 +313 4 2 7 1 21 82 98 20 +314 4 2 7 1 98 82 81 20 +315 4 2 7 1 97 98 82 68 +316 4 2 7 1 98 82 68 69 +317 4 2 7 1 68 69 82 41 +318 4 2 7 1 98 82 69 81 +319 4 2 7 1 82 41 69 81 +320 4 2 7 1 69 41 40 81 +321 4 2 7 1 16 17 60 78 +322 4 2 7 1 17 60 78 79 +323 4 2 7 1 78 79 60 98 +324 4 2 7 1 17 60 79 61 +325 4 2 7 1 60 98 79 61 +326 4 2 7 1 79 98 99 61 +327 4 2 7 1 78 79 98 36 +328 4 2 7 1 79 98 36 37 +329 4 2 7 1 36 37 98 69 +330 4 2 7 1 79 98 37 99 +331 4 2 7 1 98 69 37 99 +332 4 2 7 1 37 69 70 99 +333 4 2 7 1 60 61 20 98 +334 4 2 7 1 61 20 98 99 +335 4 2 7 1 98 99 20 81 +336 4 2 7 1 61 20 99 19 +337 4 2 7 1 20 81 99 19 +338 4 2 7 1 99 81 80 19 +339 4 2 7 1 98 99 81 69 +340 4 2 7 1 99 81 69 70 +341 4 2 7 1 69 70 81 40 +342 4 2 7 1 99 81 70 80 +343 4 2 7 1 81 40 70 80 +344 4 2 7 1 70 40 39 80 +345 4 2 7 1 17 2 61 79 +346 4 2 7 1 2 61 79 50 +347 4 2 7 1 79 50 61 99 +348 4 2 7 1 2 61 50 18 +349 4 2 7 1 61 99 50 18 +350 4 2 7 1 50 99 89 18 +351 4 2 7 1 79 50 99 37 +352 4 2 7 1 50 99 37 6 +353 4 2 7 1 37 6 99 70 +354 4 2 7 1 50 99 6 89 +355 4 2 7 1 99 70 6 89 +356 4 2 7 1 6 70 38 89 +357 4 2 7 1 61 18 19 99 +358 4 2 7 1 18 19 99 89 +359 4 2 7 1 99 89 19 80 +360 4 2 7 1 18 19 89 3 +361 4 2 7 1 19 80 89 3 +362 4 2 7 1 89 80 51 3 +363 4 2 7 1 99 89 80 70 +364 4 2 7 1 89 80 70 38 +365 4 2 7 1 70 38 80 39 +366 4 2 7 1 89 80 38 51 +367 4 2 7 1 80 39 38 51 +368 4 2 7 1 38 39 7 51 +369 2 2 1 1 1 28 49 +370 2 2 1 1 49 90 28 +371 2 2 1 1 49 90 5 +372 2 2 1 1 48 90 5 +373 2 2 1 1 90 28 4 +374 2 2 1 1 90 4 52 +375 2 2 1 1 48 90 52 +376 2 2 1 1 48 8 52 +377 2 2 4 1 27 4 52 +378 2 2 4 1 88 27 52 +379 2 2 4 1 8 88 52 +380 2 2 4 1 8 88 47 +381 2 2 4 1 88 26 27 +382 2 2 4 1 88 26 87 +383 2 2 4 1 88 87 47 +384 2 2 4 1 87 46 47 +385 2 2 4 1 25 26 87 +386 2 2 4 1 25 86 87 +387 2 2 4 1 86 46 87 +388 2 2 4 1 86 45 46 +389 2 2 4 1 24 25 86 +390 2 2 4 1 24 85 86 +391 2 2 4 1 85 45 86 +392 2 2 4 1 85 44 45 +393 2 2 4 1 24 85 23 +394 2 2 4 1 84 85 23 +395 2 2 4 1 44 85 84 +396 2 2 4 1 43 44 84 +397 2 2 4 1 84 22 23 +398 2 2 4 1 83 84 22 +399 2 2 4 1 83 43 84 +400 2 2 4 1 83 42 43 +401 2 2 4 1 83 21 22 +402 2 2 4 1 82 83 21 +403 2 2 4 1 42 83 82 +404 2 2 4 1 41 42 82 +405 2 2 4 1 82 20 21 +406 2 2 4 1 81 82 20 +407 2 2 4 1 81 41 82 +408 2 2 4 1 40 41 81 +409 2 2 4 1 81 19 20 +410 2 2 4 1 80 81 19 +411 2 2 4 1 40 81 80 +412 2 2 4 1 40 80 39 +413 2 2 4 1 80 3 19 +414 2 2 4 1 80 3 51 +415 2 2 4 1 80 51 39 +416 2 2 4 1 51 7 39 +$EndElements diff --git a/examples/Freefem/validation/3D/3D_pointFlex/comparaison_script3d_3pt.py b/examples/Freefem/validation/3D/3D_pointFlex/comparaison_script3d_3pt.py new file mode 100644 index 00000000..180807f6 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_pointFlex/comparaison_script3d_3pt.py @@ -0,0 +1,189 @@ +import json +import os +import sys +import numpy as np +import matplotlib.pyplot as plt + +from sofa_beam3d_threept import sofaRun, L, H, W +from pyfreefem import FreeFemRunner + + +def _rms(a, b): + """Discrete RMS, normalized by sqrt(n).""" + return np.linalg.norm(a - b) / np.sqrt(a.size) + + +def _default_params_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_beam3d_threept.json") + + +def _default_mesh_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "beam3d_tet.msh") + + +def _to_freefem_path(path): + """FreeFEM string literals on Windows can choke on backslashes; use forward slashes.""" + return path.replace(os.sep, "/") + + +def _pair_by_coordinates(x_a, y_a, z_a, x_b, y_b, z_b, tol=1e-6, snap=1e-6): + """ + Pair nodes between two independently-loaded copies of "the same" mesh + by (rounded) coordinates. + + The raw mesh file carries sub-tolerance floating noise from mesh + generation (e.g. x=0.8999999999997362 instead of exactly 0.9), and + this noise differs slightly across nodes that are conceptually at the + same coordinate. SOFA keeps that raw noise; FreeFEM's writer appears + to clean/round it on output. Sorting on the RAW values is therefore + unsafe: a primary sort key (x) that is only "almost tied" can order + differently between the two sources, which then contaminates the + secondary/tertiary keys (y, z) and silently mispairs otherwise + identical points. Snapping to a grid well above the noise floor but + well below the real geometric spacing fixes this. + """ + x_a, y_a, z_a = map(np.asarray, (x_a, y_a, z_a)) + x_b, y_b, z_b = map(np.asarray, (x_b, y_b, z_b)) + + print(f"[diag] n_sofa={x_a.size} n_freefem={x_b.size}") + print(f"[diag] SOFA bbox: x[{x_a.min():.6f},{x_a.max():.6f}] " + f"y[{y_a.min():.6f},{y_a.max():.6f}] z[{z_a.min():.6f},{z_a.max():.6f}]") + print(f"[diag] FreeFEM bbox: x[{x_b.min():.6f},{x_b.max():.6f}] " + f"y[{y_b.min():.6f},{y_b.max():.6f}] z[{z_b.min():.6f},{z_b.max():.6f}]") + + if x_a.size != x_b.size: + raise ValueError( + f"Node COUNT mismatch: SOFA has {x_a.size} nodes, " + f"FreeFEM has {x_b.size} nodes. Check that both loaded the " + f"exact same mesh file (same path, no stale results.txt)." + ) + + # Snap to a grid well above float noise (~1e-13) but well below the + # real mesh spacing, so ties sort consistently on both sides. + def snap_(v): + return np.round(v / snap) * snap + + xs_a, ys_a, zs_a = snap_(x_a), snap_(y_a), snap_(z_a) + xs_b, ys_b, zs_b = snap_(x_b), snap_(y_b), snap_(z_b) + + order_a = np.lexsort((zs_a, ys_a, xs_a)) + order_b = np.lexsort((zs_b, ys_b, xs_b)) + + da = np.stack([x_a[order_a], y_a[order_a], z_a[order_a]], axis=1) + db = np.stack([x_b[order_b], y_b[order_b], z_b[order_b]], axis=1) + diff = np.linalg.norm(da - db, axis=1) + worst = np.argsort(diff)[::-1][:10] + + print(f"[diag] max sorted-coordinate discrepancy = {diff.max():.6e} (tol={tol:.1e})") + if diff.max() >= tol: + print("[diag] worst offenders (sorted rank, SOFA xyz, FreeFEM xyz, dist):") + for i in worst: + print(f" rank={i:4d} sofa={tuple(np.round(da[i], 8))} " + f"ff={tuple(np.round(db[i], 8))} dist={diff[i]:.3e}") + + if not (np.allclose(x_a[order_a], x_b[order_b], atol=tol) + and np.allclose(y_a[order_a], y_b[order_b], atol=tol) + and np.allclose(z_a[order_a], z_b[order_b], atol=tol)): + raise ValueError("Node coordinates don't match between SOFA and FreeFEM meshes.") + perm = np.empty_like(order_b) + perm[order_b] = order_a + return perm + + +if __name__ == "__main__": + + config_file = sys.argv[1] if len(sys.argv) > 1 else _default_params_path() + with open(config_file) as f: + cfg = json.load(f) + + P = float(cfg["P"]) + young_modulus = float(cfg["youngModulus"]) + poisson_ratio = float(cfg["poissonRatio"]) + mesh_file = _default_mesh_path() + + os.makedirs("results", exist_ok=True) + ff_out_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), + "results", "freefem_beam3d_threept_raw.txt") + + # ====== Run FreeFEM (writes plain-text results to ff_out_path) ===== + runner = FreeFemRunner("freefem_beam3d_threept.edp") + runner.execute({ + 'P': P, + 'youngModulus': young_modulus, + 'poissonRatio': poisson_ratio, + 'meshFile': _to_freefem_path(mesh_file), + 'outFile': _to_freefem_path(ff_out_path), + }) + + if not os.path.isfile(ff_out_path): + raise RuntimeError( + f"FreeFEM did not produce the expected output file: {ff_out_path}\n" + f"Check the FreeFEM console output above for compile/runtime errors." + ) + + raw = np.loadtxt(ff_out_path) + x_ff, y_ff, z_ff = raw[:, 0], raw[:, 1], raw[:, 2] + ux_ff, uy_ff, uz_ff = raw[:, 3], raw[:, 4], raw[:, 5] + + # ========== Run SOFA =========== + pos0_sofa, u_sofa = sofaRun(mesh_file=mesh_file, P=P, + young_modulus=young_modulus, + poisson_ratio=poisson_ratio) + x_sofa, y_sofa, z_sofa = pos0_sofa[:, 0], pos0_sofa[:, 1], pos0_sofa[:, 2] + ux_sofa, uy_sofa, uz_sofa = u_sofa[:, 0], u_sofa[:, 1], u_sofa[:, 2] + + perm = _pair_by_coordinates(x_sofa, y_sofa, z_sofa, x_ff, y_ff, z_ff) + ux_ff_p = ux_ff[perm] + uy_ff_p = uy_ff[perm] + uz_ff_p = uz_ff[perm] + + rms_ux = _rms(ux_sofa, ux_ff_p) + rms_uy = _rms(uy_sofa, uy_ff_p) + rms_uz = _rms(uz_sofa, uz_ff_p) + + # Midspan deflection: node on the loaded top edge closest to (L/2, H, W/2) + mid_idx = np.argmin(np.abs(x_sofa - L / 2) + np.abs(y_sofa - H) + np.abs(z_sofa - W / 2)) + w_sofa = uy_sofa[mid_idx] + w_ff = uy_ff_p[mid_idx] + + with open("results/comparison_beam3d_threept_results.txt", 'w') as f: + header = (f"{'x':>10} {'y':>10} {'z':>10} {'ux_sofa':>12} {'ux_ff':>12} " + f"{'uy_sofa':>12} {'uy_ff':>12} {'uz_sofa':>12} {'uz_ff':>12}") + f.write(header + "\n") + f.write("-" * len(header) + "\n") + for x, y, z, uxs, uxf, uys, uyf, uzs, uzf in zip( + x_sofa, y_sofa, z_sofa, ux_sofa, ux_ff_p, uy_sofa, uy_ff_p, uz_sofa, uz_ff_p): + f.write(f"{x:10.4f} {y:10.4f} {z:10.4f} {uxs:12.6e} {uxf:12.6e} " + f"{uys:12.6e} {uyf:12.6e} {uzs:12.6e} {uzf:12.6e}\n") + + f.write("\n") + f.write("RMS norms (SOFA vs FreeFEM, discrete, sqrt(n)-normalized)\n") + f.write("-" * 55 + "\n") + f.write(f" RMS_ux (SOFA vs FF) = {rms_ux:.6e}\n") + f.write(f" RMS_uy (SOFA vs FF) = {rms_uy:.6e}\n") + f.write(f" RMS_uz (SOFA vs FF) = {rms_uz:.6e}\n") + f.write("\n") + f.write(f" w_sofa (near midspan) = {w_sofa:.6e}\n") + f.write(f" w_ff (near midspan) = {w_ff:.6e}\n") + + print(f"RMS_ux (SOFA vs FF) = {rms_ux:.6e}") + print(f"RMS_uy (SOFA vs FF) = {rms_uy:.6e}") + print(f"RMS_uz (SOFA vs FF) = {rms_uz:.6e}") + print(f"w_sofa = {w_sofa:.6e} | w_ff = {w_ff:.6e}") + + # ---- Visualization: parity plots (u_sofa vs u_ff), one per component ---- + fig, axes = plt.subplots(1, 3, figsize=(15, 5)) + for ax, (u_s, u_f, label) in zip(axes, [ + (ux_sofa, ux_ff_p, 'ux'), (uy_sofa, uy_ff_p, 'uy'), (uz_sofa, uz_ff_p, 'uz')]): + ax.scatter(u_s, u_f, s=15, alpha=0.8) + lo = min(u_s.min(), u_f.min()) + hi = max(u_s.max(), u_f.max()) + ax.plot([lo, hi], [lo, hi], 'r--', linewidth=1) + ax.set_xlabel(f'{label}_sofa') + ax.set_ylabel(f'{label}_ff') + ax.set_title(label) + + fig.suptitle("3D Beam — 3-Point Bending — SOFA vs FreeFEM (parity)", fontsize=14) + plt.tight_layout() + fig.savefig("results/comparison_beam3d_threept_fields.png", dpi=150) + plt.close(fig) \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_pointFlex/freefem_beam3d_threept.edp b/examples/Freefem/validation/3D/3D_pointFlex/freefem_beam3d_threept.edp new file mode 100644 index 00000000..fe7b233b --- /dev/null +++ b/examples/Freefem/validation/3D/3D_pointFlex/freefem_beam3d_threept.edp @@ -0,0 +1,120 @@ +load "gmsh" +load "msh3" + +DEFAULT (P, 10.0) +DEFAULT (youngModulus, 1000.0) +DEFAULT (poissonRatio, 0.3) +DEFAULT (meshFile, "beam3d_tet.msh") +DEFAULT (outFile, "freefem_beam3d_threept_out.txt") + +real Pload = $P; +real E = $youngModulus; +real nu = $poissonRatio; + +real mu = E / (2.*(1.+nu)); +real lambda = E*nu / ((1.+nu)*(1.-2.*nu)); + +mesh3 Th = gmshload3("$meshFile"); + +real Lbeam = 1.0; +real Hbeam = 0.2; +real tol = 1e-6; + +// =============== Locate support / load line nodes =============== +// (Left: x=0,y=0 - pin | Right: x=Lbeam,y=0 - roller | Load: x=Lbeam/2,y=Hbeam) +int maxN = 50; +int[int] leftIdx(maxN), rightIdx(maxN), loadIdx(maxN); +int nLeft = 0, nRight = 0, nLoad = 0; + +for (int i = 0; i < Th.nv; i++) { + real xi = Th(i).x, yi = Th(i).y; + if (abs(xi - 0.0) < tol && abs(yi - 0.0) < tol) { leftIdx[nLeft] = i; nLeft++; } + if (abs(xi - Lbeam) < tol && abs(yi - 0.0) < tol) { rightIdx[nRight] = i; nRight++; } + if (abs(xi - Lbeam/2.) < tol && abs(yi - Hbeam) < tol) { loadIdx[nLoad] = i; nLoad++; } +} +if (nLeft == 0 || nRight == 0 || nLoad == 0) { + cout << "ERROR: could not locate BC lines (nLeft=" << nLeft + << ", nRight=" << nRight << ", nLoad=" << nLoad << ")" << endl; +} + +// =============== Trapezoidal weights along z for the midspan load line =============== +real[int] zLoad(nLoad), wLoad(nLoad); +for (int k = 0; k < nLoad; k++) zLoad[k] = Th(loadIdx[k]).z; + +int[int] order(nLoad); +for (int k = 0; k < nLoad; k++) order[k] = k; +for (int a = 1; a < nLoad; a++) { + int key = order[a]; + real keyZ = zLoad[key]; + int b = a - 1; + while (b >= 0 && zLoad[order[b]] > keyZ) { + order[b+1] = order[b]; + b--; + } + order[b+1] = key; +} + +real wsum = 0.0; +for (int k = 0; k < nLoad; k++) { + real left = (k > 0) ? zLoad[order[k]] - zLoad[order[k-1]] : 0.0; + real right = (k < nLoad - 1) ? zLoad[order[k+1]] - zLoad[order[k]] : 0.0; + wLoad[order[k]] = 0.5*(left+right); + wsum += wLoad[order[k]]; +} +for (int k = 0; k < nLoad; k++) wLoad[k] = wLoad[k] / wsum; + +// =============== Vector P1 space (dofs interleaved 3*i, 3*i+1, 3*i+2 per node) =============== +fespace Wh(Th, [P1, P1, P1]); +Wh [ux, uy, uz], [vx, vy, vz]; + +varf vElasticity([ux, uy, uz], [vx, vy, vz]) = + int3d(Th)( + lambda*(dx(ux)+dy(uy)+dz(uz))*(dx(vx)+dy(vy)+dz(vz)) + + 2.*mu*( dx(ux)*dx(vx) + dy(uy)*dy(vy) + dz(uz)*dz(vz) ) + + mu*( dy(ux)+dx(uy) )*( dy(vx)+dx(vy) ) + + mu*( dz(ux)+dx(uz) )*( dz(vx)+dx(vz) ) + + mu*( dz(uy)+dy(uz) )*( dz(vy)+dy(vz) ) + ); + +matrix A = vElasticity(Wh, Wh); +real[int] b(Wh.ndof); +b = 0.0; + +// ---- Midspan distributed line load (-y direction), split by tributary length ---- +for (int k = 0; k < nLoad; k++) { + int i = loadIdx[k]; + b[3*i + 1] -= Pload * wLoad[k]; +} + +real tgv = 1e30; +// ---- Left support line: pin (ux=uy=uz=0) ---- +for (int k = 0; k < nLeft; k++) { + int i = leftIdx[k]; + A(3*i, 3*i) = tgv; b[3*i] = 0.0; + A(3*i+1, 3*i+1) = tgv; b[3*i+1] = 0.0; + A(3*i+2, 3*i+2) = tgv; b[3*i+2] = 0.0; +} +// ---- Right support line: roller (uy=0 only) ---- +for (int k = 0; k < nRight; k++) { + int i = rightIdx[k]; + A(3*i+1, 3*i+1) = tgv; b[3*i+1] = 0.0; +} + +set(A, solver=sparsesolver); +real[int] sol = A^-1 * b; +ux[] = sol; + +cout << "[sanity check] midspan deflection uy = " << uy(Lbeam/2., Hbeam, 0.1) + << " (expected: negative)" << endl; + +// =============== Write results as plain text (no pyfreefem export macros) =============== +{ + ofstream fout("$outFile"); + fout.precision(12); + for (int i = 0; i < Th.nv; i++) { + fout << Th(i).x << " " << Th(i).y << " " << Th(i).z << " " + << sol[3*i] << " " << sol[3*i+1] << " " << sol[3*i+2] << endl; + } +} + +cout << "Results written to $outFile" << endl; diff --git a/examples/Freefem/validation/3D/3D_pointFlex/params_beam3d_threept.json b/examples/Freefem/validation/3D/3D_pointFlex/params_beam3d_threept.json new file mode 100644 index 00000000..9bc56483 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_pointFlex/params_beam3d_threept.json @@ -0,0 +1,5 @@ +{ + "P": 10, + "youngModulus": 1000, + "poissonRatio": 0.3 +} \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_pointFlex/sofa_beam3d_threept.py b/examples/Freefem/validation/3D/3D_pointFlex/sofa_beam3d_threept.py new file mode 100644 index 00000000..9d68848c --- /dev/null +++ b/examples/Freefem/validation/3D/3D_pointFlex/sofa_beam3d_threept.py @@ -0,0 +1,192 @@ +""" +3D Beam Simulation - Three-Point Bending - P1 Tetrahedra +Cross-validation SOFA vs FreeFEM + +Boundary conditions (all identified by exact node-coordinate matching, +since this mesh happens to have exact nodes on x=0, x=L, and x=L/2): + - Left support line (x=0, y=0): pin -> ux = uy = uz = 0 + - Right support line (x=L, y=0): roller -> uy = 0 only + - Midspan load line (x=L/2, y=H): distributed line load along -y, + split across the nodes on that line using trapezoidal (tributary + length) weights so that the total applied force equals P. +""" +import json +import os +import sys +import numpy as np +import Sofa +import Sofa.Core +import Sofa.Simulation + +RESULTS_DIR = "results" + +L = 1.0 +H = 0.2 +W = 0.2 +TOL = 1e-6 + + +def _line_load_weights(z_values): + """Trapezoidal tributary-length weights along a 1D line (weights sum to 1).""" + z_values = np.asarray(z_values) + order = np.argsort(z_values) + z_sorted = z_values[order] + n = len(z_sorted) + w_sorted = np.zeros(n) + for k in range(n): + left = z_sorted[k] - z_sorted[k - 1] if k > 0 else 0.0 + right = z_sorted[k + 1] - z_sorted[k] if k < n - 1 else 0.0 + w_sorted[k] = 0.5 * (left + right) + w_sorted /= w_sorted.sum() + w = np.empty(n) + w[order] = w_sorted + return w + + +def _default_mesh_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "beam3d_tet.msh") + + +def _default_params_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_beam3d_threept.json") + + +def create_scene_args(rootNode, mesh_file, P, young_modulus, poisson_ratio, tol=TOL): + requiredPlugins = [ + "Elasticity", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.IO.Mesh", + "Sofa.Component.LinearSolver.Direct", + "Sofa.Component.MechanicalLoad", + "Sofa.Component.ODESolver.Backward", + "Sofa.Component.StateContainer", + "Sofa.Component.Topology.Container.Dynamic", + "Sofa.Component.Visual", + "Sofa.GL.Component.Rendering3D", + ] + rootNode.addObject('RequiredPlugin', pluginName=requiredPlugins) + rootNode.addObject('DefaultAnimationLoop') + rootNode.addObject('VisualStyle', displayFlags=["showBehaviorModels", "showForceFields"]) + + template = "Vec3d" + + with rootNode.addChild('Beam') as Beam: + Beam.addObject('NewtonRaphsonSolver' + , name="newtonSolver" + , printLog=True + , maxNbIterationsNewton=30 + , absoluteResidualStoppingThreshold=1e-12) + Beam.addObject('SparseLDLSolver' + , name="linearSolver" + , template="CompressedRowSparseMatrixd") + Beam.addObject('StaticSolver' + , name="staticSolver" + , newtonSolver="@newtonSolver" + , linearSolver="@linearSolver") + + loader = Beam.addObject('MeshGmshLoader', name="loader", filename=mesh_file) + + nodes = np.array(loader.position.value) + N = len(nodes) + x, y, z = nodes[:, 0], nodes[:, 1], nodes[:, 2] + + # ---- Left support line (x=0, y=0): pin (ux=uy=uz=0) ---- + left_idx = np.where(np.isclose(x, 0.0, atol=tol) & np.isclose(y, 0.0, atol=tol))[0] + + # ---- Right support line (x=L, y=0): roller (uy=0 only) ---- + right_idx = np.where(np.isclose(x, L, atol=tol) & np.isclose(y, 0.0, atol=tol))[0] + + # ---- Midspan load line (x=L/2, y=H): distributed line load ---- + load_idx = np.where(np.isclose(x, L / 2.0, atol=tol) & np.isclose(y, H, atol=tol))[0] + + if len(left_idx) == 0 or len(right_idx) == 0 or len(load_idx) == 0: + raise RuntimeError( + f"BC node search failed: left={len(left_idx)}, " + f"right={len(right_idx)}, load={len(load_idx)}. " + f"Check that this mesh really has exact nodes at x=0, x=L, x=L/2." + ) + + weights = _line_load_weights(z[load_idx]) + forces_y = -P * weights # applied along -y + + dofs = Beam.addObject('MechanicalObject' + , name="dofs" + , template=template + , position="@loader.position" + , showObject=True + , showObjectScale=0.01) + + Beam.addObject('TetrahedronSetTopologyContainer' + , name="topology" + , src="@loader") + Beam.addObject('TetrahedronSetTopologyModifier') + + Beam.addObject('LinearSmallStrainFEMForceField' + , name="FEM" + , template=template + , youngModulus=young_modulus + , poissonRatio=poisson_ratio + , topology="@topology") + + Beam.addObject('FixedProjectiveConstraint' + , name="leftSupport" + , indices=left_idx.tolist()) + + Beam.addObject('PartialFixedProjectiveConstraint' + , name="rightSupport" + , fixedDirections=[0, 1, 0] + , indices=right_idx.tolist()) + + Beam.addObject('ConstantForceField' + , name="MidspanLoad" + , indices=load_idx.tolist() + , forces=[[0.0, fy, 0.0] for fy in forces_y]) + + return rootNode, dofs, nodes.copy() + + +def createScene(rootNode): + with open(_default_params_path()) as f: + cfg = json.load(f) + create_scene_args(rootNode + , mesh_file=_default_mesh_path() + , P=float(cfg["P"]) + , young_modulus=float(cfg["youngModulus"]) + , poisson_ratio=float(cfg["poissonRatio"])) + return rootNode + + +def sofaRun(mesh_file, P, young_modulus, poisson_ratio): + root = Sofa.Core.Node("root") + _, dofs, pos0 = create_scene_args(root + , mesh_file=mesh_file + , P=P + , young_modulus=young_modulus + , poisson_ratio=poisson_ratio) + Sofa.Simulation.init(root) + Sofa.Simulation.animate(root, root.dt.value) + + pos_final = np.array(dofs.position.toList()) + u = pos_final[:, :3] - pos0[:, :3] + x0 = pos0[:, :3] + + os.makedirs(RESULTS_DIR, exist_ok=True) + out_path = os.path.join(RESULTS_DIR, "sofa_beam3d_threept_results.txt") + with open(out_path, 'w') as f: + f.write(f"{'x0':>12} {'y0':>12} {'z0':>12} {'ux':>12} {'uy':>12} {'uz':>12}\n") + f.write("-" * 78 + "\n") + for (xi, yi, zi), (uxi, uyi, uzi) in zip(x0, u): + f.write(f"{xi:12.6f} {yi:12.6f} {zi:12.6f} {uxi:12.6f} {uyi:12.6f} {uzi:12.6f}\n") + + return x0, u + + +if __name__ == "__main__": + config_file = sys.argv[1] if len(sys.argv) > 1 else _default_params_path() + with open(config_file) as f: + cfg = json.load(f) + + sofaRun(mesh_file=_default_mesh_path() + , P=float(cfg["P"]) + , young_modulus=float(cfg["youngModulus"]) + , poisson_ratio=float(cfg["poissonRatio"])) \ No newline at end of file