From f5a5ab225eb650e4c13b214b94cce2b95bcf5b80 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Thu, 10 Sep 2026 09:08:02 +1000 Subject: [PATCH 1/4] Added new reffe tutorials --- models/clamped_disk.msh | 2635 ++++++++++++++++++++++++++++++ models/lshape.msh | 1680 +++++++++++++++++++ src/argyris_biharmonic.jl | 125 ++ src/arnold_winther_elasticity.jl | 158 ++ src/gls_stokes.jl | 130 ++ src/hhj_plate.jl | 166 ++ src/morley_biharmonic.jl | 128 ++ src/mtw_darcy_stokes.jl | 126 ++ src/regge_metric.jl | 122 ++ 9 files changed, 5270 insertions(+) create mode 100644 models/clamped_disk.msh create mode 100644 models/lshape.msh create mode 100644 src/argyris_biharmonic.jl create mode 100644 src/arnold_winther_elasticity.jl create mode 100644 src/gls_stokes.jl create mode 100644 src/hhj_plate.jl create mode 100644 src/morley_biharmonic.jl create mode 100644 src/mtw_darcy_stokes.jl create mode 100644 src/regge_metric.jl diff --git a/models/clamped_disk.msh b/models/clamped_disk.msh new file mode 100644 index 0000000..57403f0 --- /dev/null +++ b/models/clamped_disk.msh @@ -0,0 +1,2635 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +2 +1 1 "boundary" +2 2 "domain" +$EndPhysicalNames +$Entities +5 4 1 0 +1 0 0 0 0 +2 1 0 0 0 +3 0 1 0 0 +4 -1 0 0 0 +5 0 -1 0 0 +1 5.551115123125783e-17 0 0 1 1 0 1 1 2 2 -3 +2 -1 5.551115123125783e-17 0 0 1 0 1 1 2 3 -4 +3 -1 -1 0 -5.551115123125783e-17 0 0 1 1 2 4 -5 +4 0 -1 0 1 -5.551115123125783e-17 0 1 1 2 5 -2 +1 -1 -1 0 1 1 0 1 2 4 1 2 3 4 +$EndEntities +$Nodes +9 649 1 649 +0 2 0 1 +1 +1 0 0 +0 3 0 1 +2 +0 1 0 +0 4 0 1 +3 +-1 0 0 +0 5 0 1 +4 +0 -1 0 +1 1 0 19 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +0.9969173337185187 0.07845909591347422 0 +0.9876883405347197 0.1564344654216953 0 +0.9723699202529145 0.2334453644588833 0 +0.9510565160292614 0.3090169951932792 0 +0.9238795320809745 0.3826834334039555 0 +0.8910065235789868 0.4539905009355246 0 +0.8526401635399679 0.522498566044479 0 +0.8090169933233242 0.5877852537399083 0 +0.7604059642573492 0.649448049902262 0 +0.707106779570334 0.7071067828027611 0 +0.6494480467864843 0.7604059669184748 0 +0.5877852508817173 0.8090169953999216 0 +0.5224985633033551 0.8526401652197309 0 +0.4539904984339875 0.8910065248535836 0 +0.3826834312304585 0.9238795329812663 0 +0.3090169934436383 0.9510565165977543 0 +0.2334453631624346 0.9723699205641642 0 +0.1564344645361618 0.9876883406749745 0 +0.07845909547534838 0.9969173337529998 0 +1 2 0 19 +24 +25 +26 +27 +28 +29 +30 +31 +32 +33 +34 +35 +36 +37 +38 +39 +40 +41 +42 +-0.07845909591347422 0.9969173337185187 0 +-0.1564344654216953 0.9876883405347197 0 +-0.2334453644588833 0.9723699202529145 0 +-0.3090169951932792 0.9510565160292614 0 +-0.3826834334039555 0.9238795320809745 0 +-0.4539905009355246 0.8910065235789868 0 +-0.522498566044479 0.8526401635399679 0 +-0.5877852537399083 0.8090169933233242 0 +-0.649448049902262 0.7604059642573492 0 +-0.7071067828027611 0.707106779570334 0 +-0.7604059669184748 0.6494480467864843 0 +-0.8090169953999216 0.5877852508817173 0 +-0.8526401652197309 0.5224985633033551 0 +-0.8910065248535836 0.4539904984339875 0 +-0.9238795329812663 0.3826834312304585 0 +-0.9510565165977543 0.3090169934436383 0 +-0.9723699205641642 0.2334453631624346 0 +-0.9876883406749745 0.1564344645361618 0 +-0.9969173337529998 0.07845909547534838 0 +1 3 0 19 +43 +44 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +-0.9969173337185187 -0.07845909591347422 0 +-0.9876883405347197 -0.1564344654216953 0 +-0.9723699202529145 -0.2334453644588833 0 +-0.9510565160292614 -0.3090169951932792 0 +-0.9238795320809745 -0.3826834334039555 0 +-0.8910065235789868 -0.4539905009355246 0 +-0.8526401635399679 -0.522498566044479 0 +-0.8090169933233242 -0.5877852537399083 0 +-0.7604059642573492 -0.649448049902262 0 +-0.707106779570334 -0.7071067828027611 0 +-0.6494480467864843 -0.7604059669184748 0 +-0.5877852508817173 -0.8090169953999216 0 +-0.5224985633033551 -0.8526401652197309 0 +-0.4539904984339875 -0.8910065248535836 0 +-0.3826834312304585 -0.9238795329812663 0 +-0.3090169934436383 -0.9510565165977543 0 +-0.2334453631624346 -0.9723699205641642 0 +-0.1564344645361618 -0.9876883406749745 0 +-0.07845909547534838 -0.9969173337529998 0 +1 4 0 19 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +0.07845909591347422 -0.9969173337185187 0 +0.1564344654216953 -0.9876883405347197 0 +0.2334453644588833 -0.9723699202529145 0 +0.3090169951932792 -0.9510565160292614 0 +0.3826834334039555 -0.9238795320809745 0 +0.4539905009355246 -0.8910065235789868 0 +0.522498566044479 -0.8526401635399679 0 +0.5877852537399083 -0.8090169933233242 0 +0.649448049902262 -0.7604059642573492 0 +0.7071067828027611 -0.707106779570334 0 +0.7604059669184748 -0.6494480467864843 0 +0.8090169953999216 -0.5877852508817173 0 +0.8526401652197309 -0.5224985633033551 0 +0.8910065248535836 -0.4539904984339875 0 +0.9238795329812663 -0.3826834312304585 0 +0.9510565165977543 -0.3090169934436383 0 +0.9723699205641642 -0.2334453631624346 0 +0.9876883406749745 -0.1564344645361618 0 +0.9969173337529998 -0.07845909547534838 0 +2 1 0 569 +81 +82 +83 +84 +85 +86 +87 +88 +89 +90 +91 +92 +93 +94 +95 +96 +97 +98 +99 +100 +101 +102 +103 +104 +105 +106 +107 +108 +109 +110 +111 +112 +113 +114 +115 +116 +117 +118 +119 +120 +121 +122 +123 +124 +125 +126 +127 +128 +129 +130 +131 +132 +133 +134 +135 +136 +137 +138 +139 +140 +141 +142 +143 +144 +145 +146 +147 +148 +149 +150 +151 +152 +153 +154 +155 +156 +157 +158 +159 +160 +161 +162 +163 +164 +165 +166 +167 +168 +169 +170 +171 +172 +173 +174 +175 +176 +177 +178 +179 +180 +181 +182 +183 +184 +185 +186 +187 +188 +189 +190 +191 +192 +193 +194 +195 +196 +197 +198 +199 +200 +201 +202 +203 +204 +205 +206 +207 +208 +209 +210 +211 +212 +213 +214 +215 +216 +217 +218 +219 +220 +221 +222 +223 +224 +225 +226 +227 +228 +229 +230 +231 +232 +233 +234 +235 +236 +237 +238 +239 +240 +241 +242 +243 +244 +245 +246 +247 +248 +249 +250 +251 +252 +253 +254 +255 +256 +257 +258 +259 +260 +261 +262 +263 +264 +265 +266 +267 +268 +269 +270 +271 +272 +273 +274 +275 +276 +277 +278 +279 +280 +281 +282 +283 +284 +285 +286 +287 +288 +289 +290 +291 +292 +293 +294 +295 +296 +297 +298 +299 +300 +301 +302 +303 +304 +305 +306 +307 +308 +309 +310 +311 +312 +313 +314 +315 +316 +317 +318 +319 +320 +321 +322 +323 +324 +325 +326 +327 +328 +329 +330 +331 +332 +333 +334 +335 +336 +337 +338 +339 +340 +341 +342 +343 +344 +345 +346 +347 +348 +349 +350 +351 +352 +353 +354 +355 +356 +357 +358 +359 +360 +361 +362 +363 +364 +365 +366 +367 +368 +369 +370 +371 +372 +373 +374 +375 +376 +377 +378 +379 +380 +381 +382 +383 +384 +385 +386 +387 +388 +389 +390 +391 +392 +393 +394 +395 +396 +397 +398 +399 +400 +401 +402 +403 +404 +405 +406 +407 +408 +409 +410 +411 +412 +413 +414 +415 +416 +417 +418 +419 +420 +421 +422 +423 +424 +425 +426 +427 +428 +429 +430 +431 +432 +433 +434 +435 +436 +437 +438 +439 +440 +441 +442 +443 +444 +445 +446 +447 +448 +449 +450 +451 +452 +453 +454 +455 +456 +457 +458 +459 +460 +461 +462 +463 +464 +465 +466 +467 +468 +469 +470 +471 +472 +473 +474 +475 +476 +477 +478 +479 +480 +481 +482 +483 +484 +485 +486 +487 +488 +489 +490 +491 +492 +493 +494 +495 +496 +497 +498 +499 +500 +501 +502 +503 +504 +505 +506 +507 +508 +509 +510 +511 +512 +513 +514 +515 +516 +517 +518 +519 +520 +521 +522 +523 +524 +525 +526 +527 +528 +529 +530 +531 +532 +533 +534 +535 +536 +537 +538 +539 +540 +541 +542 +543 +544 +545 +546 +547 +548 +549 +550 +551 +552 +553 +554 +555 +556 +557 +558 +559 +560 +561 +562 +563 +564 +565 +566 +567 +568 +569 +570 +571 +572 +573 +574 +575 +576 +577 +578 +579 +580 +581 +582 +583 +584 +585 +586 +587 +588 +589 +590 +591 +592 +593 +594 +595 +596 +597 +598 +599 +600 +601 +602 +603 +604 +605 +606 +607 +608 +609 +610 +611 +612 +613 +614 +615 +616 +617 +618 +619 +620 +621 +622 +623 +624 +625 +626 +627 +628 +629 +630 +631 +632 +633 +634 +635 +636 +637 +638 +639 +640 +641 +642 +643 +644 +645 +646 +647 +648 +649 +0.6850115571317198 -0.6339723868751211 0 +-0.6838224475295513 0.6321189657422444 0 +0.6329856152397908 0.691093509089582 0 +-0.6321189657422444 -0.6838224475295513 0 +-0.9277172783897996 -0.03724384464503164 0 +-0.04266513669016786 0.9219826058616457 0 +0.03329274234547587 -0.9268538419878467 0 +0.9305110968302381 0.03655988065304996 0 +-0.3866992990735282 0.8485719001721408 0 +-0.8381155400980319 -0.3853331934609263 0 +0.3873413998849349 -0.8436040686140578 0 +0.8424928462114599 0.3897649285727175 0 +0.3178655140766905 0.8807446615490984 0 +-0.8730947568819056 0.3267811291476899 0 +-0.3277512381816218 -0.8822751151187025 0 +0.8744031165281956 -0.3340974361193404 0 +-0.5748831533830876 0.7408808562115882 0 +0.7379816774327458 0.5716520798905029 0 +-0.7313099308818299 -0.5765182657746002 0 +0.5759785222048306 -0.7303427640136195 0 +-0.1806938342787924 0.9160291714110738 0 +-0.9094398970480013 -0.1788182471872823 0 +0.1847436790953905 -0.9283489568380573 0 +0.9021985985771287 0.1808340271974244 0 +-0.7825823877482552 0.514429844953525 0 +0.5154611350200825 0.7762542770954921 0 +-0.5232149375194273 -0.776483304424457 0 +-0.9115948402927155 0.1830606860414794 0 +0.1795990314301271 0.9166293790248151 0 +-0.1768703810288315 -0.9220059734504649 0 +0.9150319140516833 -0.1738752843771139 0 +0.8143084799139373 -0.4457029237159921 0 +-0.6285642397987277 0.68737545610038 0 +-0.6058914127710918 0.6144151562410898 0 +-0.6602626462414919 0.5575317684720412 0 +-0.5831286225906004 0.5412418606847709 0 +-0.6347120602747737 0.4820051316030328 0 +-0.5586778203068314 0.4667190510908324 0 +-0.5065358303990375 0.5254194912463788 0 +-0.4819944732963296 0.4508551554173488 0 +-0.5343046641309648 0.3922972364179117 0 +-0.457484258123716 0.3763182292541459 0 +-0.4051689867226068 0.4348716659906599 0 +-0.3806146442669688 0.3602980956726042 0 +-0.4329227634437864 0.3017387122448661 0 +-0.3560544911522563 0.2857196937418086 0 +-0.3037475246277082 0.3442800376486186 0 +-0.2791865298891108 0.2697009350274199 0 +-0.3314934451166035 0.2111405489356208 0 +-0.2546252525915447 0.1951215982244068 0 +-0.2023182985655865 0.2536819528010761 0 +-0.1777570287342576 0.1791026216084721 0 +-0.2300639886718802 0.1205422725959674 0 +-0.1531957602617488 0.1045232920400485 0 +-0.1008888001445425 0.163083640857133 0 +-0.07632753055414659 0.08850431044427176 0 +-0.1286344905193996 0.02994396153989573 0 +-0.05176626074485299 0.01392497987054655 0 +0.0005406992268303067 0.07248532879617614 0 +0.02510196905192237 -0.002094001797194897 0 +-0.02720499092476485 -0.06065435071741808 0 +0.04966323887991453 -0.0766733323885343 0 +0.07740892903036597 0.05646634712682115 0 +-0.002643721099341762 -0.1352336813091547 0 +0.0742245087102066 -0.1512526629826178 0 +0.02191754872903133 -0.2098130119023298 0 +-0.0549506810749555 -0.1937940302252365 0 +0.09878577853916426 -0.2258319935814954 0 +0.151092738524047 -0.1672716446571785 0 +-0.03038941125031921 -0.2683733608164492 0 +-0.1072576410475373 -0.2523543791410817 0 +0.1756540083564924 -0.2418509752634698 0 +0.2279609683438383 -0.183290626336457 0 +-0.08269637122971006 -0.3269337097234334 0 +-0.1595646010177907 -0.3109147280515776 0 +0.2525222381773556 -0.2578699569441774 0 +0.3048291981644234 -0.1993096080154608 0 +0.2002152781841685 -0.3164303058801503 0 +-0.1350033312016815 -0.3854940586218341 0 +-0.2118715609785584 -0.369475076952708 0 +-0.2364328307933255 -0.2948957463854877 0 +0.2770835080268791 -0.33244928755762 0 +0.2247765480112749 -0.3910096365126564 0 +-0.2887397907517605 -0.3534560952843729 0 +-0.3133010605697115 -0.2788767647226793 0 +-0.1873102911623076 -0.4440544075133316 0 +-0.110442061386904 -0.4600733891885752 0 +-0.3656080205260602 -0.3374371136191483 0 +-0.3901692903528933 -0.2628577830460486 0 +0.3016447778767326 -0.4070286181905596 0 +0.2493378178470395 -0.4655889671587843 0 +-0.4424762503147092 -0.3214181319542584 0 +-0.4670375201525501 -0.2468388013663505 0 +0.2802679283326912 -0.1247302774110381 0 +0.3571361581510669 -0.1407492590893701 0 +0.3325748883203158 -0.06616992848419528 0 +0.409443118144836 -0.08218891016667836 0 +0.3848818483102479 -0.007609579551478928 0 +0.4340043880239631 -0.1567682407812574 0 +0.4617500781399012 -0.02362856124027516 0 +0.4371888083029124 0.05095076938206623 0 +0.5140570381338003 0.03493178769693001 0 +0.4894957682949917 0.1095111183152835 0 +0.4126275384679907 0.1255300999934463 0 +0.4649344984580241 0.1840904489259752 0 +0.3880662686403777 0.2001094305967257 0 +0.4403732286117812 0.2586697795184731 0 +0.3635049988220636 0.2746887611830149 0 +0.4158119588041229 0.3332491101589672 0 +0.3389437290065473 0.3492680917613007 0 +-0.1627490213543382 -0.5186337380790543 0 +-0.08588079158007235 -0.5346527197582308 0 +0.3912506889547761 0.4078284406501589 0 +0.3143824591827241 0.4238474223163172 0 +0.2620754992030768 0.365287073418868 0 +0.2375142294027929 0.4398664039862253 0 +-0.4179149804877122 -0.3959974625260022 0 +0.5418027282839577 0.1680714672843636 0 +0.2898211894679154 0.4984267529996396 0 +0.2129529595832663 0.5144457345585189 0 +0.1606459996057868 0.4558853856568906 0 +0.1360847297788346 0.5304647162568178 0 +0.08377776976719387 0.4719043673309522 0 +0.05921649995744732 0.5464836979493668 0 +-0.4147305601593946 -0.1882784524324108 0 +-0.4915987899852056 -0.1722594707570037 0 +0.326206047717301 -0.4816079488435228 0 +0.3785130077395578 -0.4230475998764114 0 +0.2738990876943166 -0.5401682978162383 0 +-0.5439057499820107 -0.2308198196966391 0 +-0.5684670198284159 -0.1562404890714662 0 +-0.5161600598157402 -0.09768014011979678 0 +-0.5930282896698604 -0.08166115843414919 0 +0.4926801886300614 0.3172301286108647 0 +0.006909540078897037 0.4879233490811522 0 +-0.01765172985868166 0.5625026796407158 0 +0.1083390395775682 0.3973250367423597 0 +0.5386183079760892 -0.0396475429276044 0 +0.590925267969862 0.01891280601353714 0 +-0.540721329657766 -0.02310080948453733 0 +-0.6175895595054711 -0.007081827792389914 0 +0.6154865378114244 -0.05566652461683589 0 +0.03465523013632243 0.6210630285805842 0 +-0.04221299970803828 0.6370820102713124 0 +-0.09451995972889021 0.5785216613165722 0 +-0.1381877515468054 -0.5932130686599366 0 +-0.06131952184982718 -0.6092320504617729 0 +-0.2150559813107487 -0.57719408697313 0 +-0.4947832102660088 -0.3799784808615445 0 +0.1883916899191631 0.5890250652893602 0 +-0.1190812295677135 0.6531009919720034 0 +-0.6207739798345173 -0.2148008380247715 0 +-0.4702219404324764 -0.4545578114135029 0 +0.4030742775843437 -0.4976269305262698 0 +0.4553812375951591 -0.4390665815676206 0 +-0.171388189594941 0.5945406430162963 0 +0.3111980388179372 0.2161284122507089 0 +-0.2609941005850525 -0.2203164158049924 0 +-0.5652825994974446 0.05147852114941959 0 +-0.6421508293320576 0.06749750284175168 0 +-0.4638530998061521 -0.03911979117785895 0 +0.05284765919871784 0.1310456777280122 0 +-0.03357383162353398 -0.4760923708840552 0 +-0.6698965195223084 -0.06564217673746783 0 +0.6677934978199362 0.002893824329439715 0 +-0.1904947115306828 -0.651773417531756 0 +-0.2673629412687072 -0.635754435846154 0 +-0.1959494594345969 0.6691199736796265 0 +-0.2482564194663332 0.6105596247197895 0 +0.3080136184867933 0.008409402116994676 0 +-0.5898438693356758 0.1260578517761204 0 +-0.6667123599998175 0.1420767475528438 0 +0.1970308578079312 -0.5241493161230879 0 +0.2215921276517784 -0.5987286467978765 0 +0.2984603575461638 -0.6147476284874858 0 +-0.1468269203608791 0.5199613128260622 0 +0.6923547676519477 -0.07168550630382665 0 +0.6400478076555162 -0.1302458552350728 0 +0.01009396030704363 0.6956423592132923 0 +0.2461533975067265 -0.6733079774879618 0 +0.1692851675982512 -0.6572889957789474 0 +-0.7189143357939557 0.08355072520677592 0 +0.479942507441026 -0.5136459122071368 0 +0.5322494674402711 -0.4550855632490189 0 +0.5076881975977492 -0.3805062326301633 0 +0.5845564274269063 -0.3965252143156355 0 +-0.5470901701918099 -0.4385388297711046 0 +0.7169160374773131 -0.1462648369236695 0 +-0.6144051826393602 0.2006371680553678 0 +-0.6912736805311674 0.2166560614476117 0 +0.1297158890120432 0.1150266960517909 0 +-0.393353710660334 -0.4705767930752237 0 +-0.4456606705750698 -0.5291371419654307 0 +0.6646090774889997 -0.2048251858441903 0 +0.1051546191759698 0.1896060266535785 0 +0.1479083181725615 -0.3749906548146686 0 +-0.5129756395031551 0.1100388700926254 0 +0.7692229974769411 -0.08770448800651892 0 +-0.2728176893067543 0.6851389553847552 0 +-0.3251247748389124 0.62657862381463 0 +-0.1436424994010182 0.7276803226235957 0 +-0.2428016714783938 -0.7103337663882238 0 +-0.1659334417890399 -0.7263527480994092 0 +-0.3196699012075533 -0.6943147846986758 0 +0.6091176972638311 -0.4711045449135956 0 +0.1938464374514617 -0.7318683264814917 0 +0.1169782075405575 -0.7158493447649176 0 +0.3230216274036498 -0.6893269591685083 0 +0.6614246572392593 -0.4125441959929385 0 +0.7414773072950991 -0.2208441675330379 0 +-0.6389664707857362 0.2752164958093252 0 +-0.715834700996223 0.291235474676964 0 +-0.3005636227606501 0.5519993818337842 0 +0.5877408476692919 -0.1888062041622193 0 +0.08696219018232422 0.6796233775323818 0 +0.06240092037475673 0.7542027081773743 0 +0.1392691502843973 0.7381837265124053 0 +-0.131818910869967 -0.177775048555968 0 +0.1542771588426681 0.04044736545287182 0 +0.0282863893630416 0.2056250083295802 0 +0.2065841188255361 0.09900771437702129 0 +0.6432322279811729 0.07747315497413447 0 +0.7201004578421735 0.06145417328417795 0 +0.6955391880059915 0.1360335039516045 0 +0.4276355474366361 -0.5722062611752607 0 +-0.005828141442080313 -0.3429526913992277 0 +-0.5962127099655523 -0.2893801686397002 0 +-0.6730809398384849 -0.2733611869870146 0 +-0.64851966997571 -0.347940517619063 0 +-0.6976422097041879 -0.1987818563423482 0 +-0.7499491697306843 -0.2573422053232097 0 +-0.01446730951333985 0.7702216898542297 0 +0.09241693772062035 -0.6412700140640266 0 +0.04010997767111849 -0.6998303630423371 0 +0.5032589119547006 -0.5871133109227291 0 +0.4519893397405351 -0.6466002698489102 0 +-0.3498170610535326 0.7005287466861612 0 +-0.4020151442252096 0.6424927083854457 0 +-0.4263296229464293 0.7159119717356837 0 +0.06467124748828876 -0.77440969374095 0 +-0.01241282267923503 -0.7594264370810861 0 +0.7724074178909001 0.1200145222547673 0 +0.7478461480382608 0.1945938529374552 0 +0.6709779181549106 0.2106128346283299 0 +0.72328487818002 0.2691731836133771 0 +0.6464166482914496 0.2851921653071502 0 +0.6987236083141367 0.3437525142718644 0 +0.6218553784628615 0.3597714959643857 0 +0.6741623384567705 0.4183318449042579 0 +0.5972941086410943 0.434350826585277 0 +-0.7682299852780666 0.2326261994913391 0 +0.6848194686851665 -0.4871054000436745 0 +0.7957407975860692 0.2517010922272693 0 +-0.7732443383237382 -0.1823459085236271 0 +-0.2968322245508992 0.7583592342258677 0 +0.7951170384552732 0.04482538645672966 0 +0.03710671149348182 0.8252649783371112 0 +-0.2182404016965317 -0.7849130969412206 0 +-0.1407321853852755 -0.8006102423639977 0 +-0.4771290438910548 0.6568937607899068 0 +-0.6635276906927164 0.3497958415635248 0 +-0.7404144347386401 0.3657259816818466 0 +-0.3442311709925586 -0.6197354541784527 0 +-0.3965001820766135 -0.6782833052971046 0 +-0.3737275753410152 -0.7512444812738379 0 +0.6891703473197399 -0.2794045164512761 0 +-0.5620982043996199 0.259197534084044 0 +0.5599951575805702 -0.3219458837048285 0 +0.48312692776607 -0.3059269020091598 0 +0.08059334934758269 0.2641853572485265 0 +0.003725119583491662 0.2802043388798078 0 +0.1574615791622629 0.2481663755717471 0 +0.6460359960336169 0.4917370873429794 0 +0.5721386600470618 0.5087344757973359 0 +0.5203268490671287 0.4503371947092306 0 +0.4957490743092566 0.5249110896642337 0 +0.5475489941474618 0.5810506645501097 0 +0.471180320998407 0.5991123239518659 0 +0.418956620458424 0.5408924213220379 0 +0.394407980024738 0.6154654769068529 0 +0.4460977030223686 0.6710528239087499 0 +0.3697619237639763 0.6896039565093771 0 +0.3175437055710231 0.6314656328140045 0 +0.2929830949632776 0.7060418257158877 0 +-0.8378019219621746 -0.234661844588577 0 +0.7628099177328946 -0.2925396123577997 0 +0.3444598745369143 0.7606174443863557 0 +0.2674681815499002 0.7760421189972208 0 +0.7251014144984641 0.4872168738548309 0 +0.8124907690292769 -0.2347924108030475 0 +0.6275722967874788 -0.5416239933155561 0 +0.8597409497649626 0.0995284733070743 0 +0.01124116866721167 -0.8384927640283142 0 +0.08923251732142623 -0.8489890244382727 0 +-0.04747291049474315 -0.8946656119081228 0 +-0.03902857936312661 0.8448010204981699 0 +-0.8022561297503157 -0.3159025543432218 0 +0.5317497686346304 -0.6568975830439494 0 +0.4794844203677769 -0.7194364212918167 0 +0.6082443252627137 -0.6562100770573966 0 +-0.06579898417029362 -0.8231654117873435 0 +0.7371123165698125 -0.429016799902048 0 +-0.6869335244871453 0.4232255750737971 0 +-0.5291479196198552 0.5983030258250597 0 +-0.4297173299201418 0.5094401908099555 0 +-0.6112207188617303 0.408356192408228 0 +-0.328307504546573 0.4188583004497218 0 +-0.2514405342687017 0.402840362091778 0 +-0.3528653198054038 0.4934346410795983 0 +-0.2268795691714013 0.3282612868020039 0 +-0.1500113395943663 0.3122423031186291 0 +-0.1745725594559795 0.3868215923017605 0 +-0.09770437134610423 0.3708026448404569 0 +-0.5097832173825091 0.3177511895844529 0 +0.2406910351111836 0.6475512913959677 0 +-0.4083616297781164 0.2271594934462066 0 +-0.3838004273720357 0.1525802186277073 0 +-0.03679422549918596 -0.6839840019987162 0 +0.2652547882122922 0.5729835965343089 0 +0.1638291509279575 0.6635987087512258 0 +0.2159694538904577 0.7213602911480597 0 +0.190750781282353 0.7927840987111272 0 +0.2424237873449288 0.8515372096920518 0 +0.1138135207851916 0.8084714789401538 0 +0.09215202179007131 0.882845620026391 0 +0.6219597714033919 0.5649233700352475 0 +0.5941096882880916 0.226631816291738 0 +-0.421093075930042 -0.6037143895449624 0 +-0.4724447663987429 -0.6602120856394974 0 +-0.3687913866889432 -0.5451557764663626 0 +-0.3164853051943518 -0.4865957168816888 0 +-0.7933913576986937 0.3077886591001269 0 +-0.8114422333240623 0.384026061610883 0 +-0.8462695584750293 0.2468539597274229 0 +-0.8229082851610032 0.1717399377251061 0 +-0.8720942116522729 0.1081744813563053 0 +-0.7438826282692289 0.1577318638372725 0 +-0.5225289003396553 -0.5131181603054938 0 +-0.5993971301046231 -0.497099178681212 0 +-0.5742718363879196 -0.569537168180805 0 +-0.6581398449786955 -0.5572323260299051 0 +-0.6762653598783518 -0.4810801970565976 0 +-0.7626924962554033 -0.4948254483119621 0 +0.05603207953710278 0.3387646878248658 0 +-0.02083614861087942 0.3547836684615803 0 +0.5449706436862314 0.3757850420539147 0 +0.4680996630475159 0.3918031176918993 0 +0.5695456676280146 0.3012102410163366 0 +-0.2236951914027011 0.5359803121019685 0 +-0.2396172218335215 -0.5026147467696063 0 +0.1447238977475731 -0.5827096650922178 0 +0.1201626279462091 -0.5081303344102408 0 +0.06785566794754563 -0.5666906833273383 0 +0.789618380641816 -0.1629892910562142 0 +0.8537167295182885 -0.1134547720468258 0 +-0.2964610146326324 -0.7654631881589586 0 +-0.2641784867733712 -0.4280354145950375 0 +0.3603205784719475 0.06696975105235027 0 +0.5172409999833216 0.2426506469412919 0 +-0.06995868977214287 0.5039423311077368 0 +0.2033996985124792 -0.1087112957360611 0 +-0.488414369681214 0.03545953945239783 0 +-0.4115461398466199 0.01944055773443783 0 +-0.3869848699737709 -0.05513877286025163 0 +-0.3346779099879514 0.003421576038355773 0 +-0.3101166401431291 -0.07115775453652214 0 +-0.2578096801962322 -0.01259740565537182 0 +-0.2332484103902417 -0.08717673627249908 0 +-0.3624236001536676 -0.1297181034856874 0 +-0.4392918299823219 -0.1136991218055015 0 +-0.285555370345205 -0.1457370851662166 0 +-0.4977123066700107 -0.5869961308861147 0 +-0.5572310020007633 -0.651763377149321 0 +0.1820228489948313 0.1735870449765938 0 +0.2588910788112543 0.1575680633142005 0 +0.2343298089934586 0.232147393910679 0 +0.2097685391667538 0.3067267244997794 0 +0.3357593086407541 0.1415490816525484 0 +-0.3069322072654493 0.1365612449664266 0 +0.4308199677620399 -0.3644872509302617 0 +0.4062586979160686 -0.2899079203048601 0 +0.4585656579446122 -0.2313475713933424 0 +0.3539517378961448 -0.3484682692335135 0 +0.5631795778333698 -0.1142268735476729 0 +-0.04858184033029572 0.2216439898994853 0 +0.7129967455437654 -0.3535788031479516 0 +0.7972410075418667 -0.3633536586620921 0 +-0.079511950889742 -0.1192146996450785 0 +0.3816974279963246 -0.2153285896954814 0 +-0.694440335454654 0.008942860684316892 0 +-0.7467618403970682 -0.0496222439031183 0 +-0.104073220692482 -0.04463536905223457 0 +0.2311453886618287 0.02442838378304429 0 +0.1724695879856952 -0.4495699854599134 0 +0.09560135816181071 -0.4335510037311469 0 +0.5222265883324517 0.6455355813777738 0 +0.04647881855152668 -0.2843923424966316 0 +-0.7672462392299269 0.02636337049639564 0 +-0.8313679391191853 -0.02625180248283805 0 +0.4333864437177134 0.7503854255434206 0 +0.3904553026189821 0.8310435186068104 0 +0.1115232484836749 0.605043099059956 0 +0.1265314686921242 -0.09269231406028598 0 +0.1233470483650975 -0.3004113241881281 0 +-0.341046715732265 -0.4120164326635789 0 +0.366682792804636 0.4823951549007959 0 +0.3293904679981248 -0.2738889386197135 0 +0.486311347999163 -0.09820789186018393 0 +0.7443531194554153 -0.01322679155771207 0 +-0.02402057054006584 0.1470646593424431 0 +0.5663639981344233 0.09349213665077094 0 +0.2834523486496857 0.08298873271602655 0 +-0.2919240012148051 -0.5611750361858989 0 +-0.5193444801197722 -0.3053991502876356 0 +-0.5375369123325764 0.1846182001907297 0 +-0.1841258708217276 -0.2363353974821059 0 +0.5550463649350516 -0.5282549646206969 0 +0.6367409087723115 -0.3378973580622271 0 +0.1019701988589556 -0.01811298346785698 0 +0.3507673175694827 -0.5561872794940546 0 +-0.645335249682332 -0.1402215073854787 0 +0.8212850660055417 -0.02975593159677046 0 +0.8924664802807598 -0.03819163435522531 0 +0.777278202145461 0.3280663671037013 0 +0.3421111794231714 0.5569381725796066 0 +0.6186709581414245 0.1520524856301569 0 +0.2557066584972785 -0.05015094681224792 0 +0.2866367690016395 0.2907077428373918 0 +-0.2204196068385771 0.7434727957089813 0 +-0.2452026031355253 0.8129065049469376 0 +0.07104008835606315 -0.35897167309959 0 +0.4435242912600473 0.4663691089121643 0 +0.1788384286775557 -0.03413196514008932 0 +0.2707146673660698 -0.7478873081895187 0 +0.220846973138834 -0.8073043959732287 0 +0.2952759372491103 -0.8224666389182197 0 +-0.7219919778407919 -0.1241328728114375 0 +0.7525721424844377 0.4033053652697941 0 +0.8072068176389302 0.4634028343627432 0 +-0.3378623303609874 -0.2042974341096725 0 +-0.7244529744654216 -0.3274352295868007 0 +0.6122817044295614 -0.2633742835532286 0 +-0.05813510142607449 -0.4015130402974772 0 +-0.06677426954117451 0.7116613409088222 0 +0.142489895517452 -0.7879616478071184 0 +-0.799071709415117 -0.1081835439765849 0 +0.8247143779349086 0.1785748712348603 0 +-0.5716514400534947 -0.363959499212938 0 +0.1852072693858891 0.3813060550746823 0 +0.3501943440265975 -0.7643823184951435 0 +0.3752940078322807 -0.6307357231642702 0 +-0.6239583999663201 -0.4225198481573567 0 +-0.7048839881143872 -0.4071106470976995 0 +0.1329003093611437 0.3227457061523159 0 +0.5354304855478919 -0.2473646778828767 0 +-0.1689658792266443 0.7977010833242698 0 +0.04329439816509172 -0.4921113526419166 0 +0.510872050836382 -0.1727869099379254 0 +0.03147080966249443 0.4133440184055105 0 +-0.09146255770264723 0.7854809098887294 0 +-0.1254500695540248 0.2376629714327017 0 +0.0187331283652128 -0.4175320220089022 0 +-0.4852285678315706 0.2431773914219676 0 +-0.2086871406124836 -0.1617560669177625 0 +-0.1563801806437185 -0.1031957179910133 0 +0.4008177190549073 -0.7046907812063005 0 +-0.1809414504513376 -0.02861638736317604 0 +-0.009012561832253432 -0.5506717015632967 0 +-0.2055027200182343 0.04596294300385304 0 +0.01554271235966964 -0.6252798024095811 0 +-0.1136324774102765 -0.6678211694344746 0 +-0.2823709476619286 0.06198192291638902 0 +-0.1178564870049736 0.8636606198978884 0 +-0.3592391736590826 0.07800090144202315 0 +-0.4606684311063275 0.1685990103580098 0 +-0.4361073635280724 0.09401984961786687 0 +-0.08921735615555107 -0.7435600017941711 0 +-0.07314310830869647 0.2962233194387768 0 +-0.04539741826731289 0.4293629991951549 0 +-0.276001009416255 0.4774189936975784 0 +-0.1222656404859865 0.4453819741499652 0 +-0.1991336425650838 0.4614007578615188 0 +-0.5866501610422127 0.3337690816445811 0 +-0.4541352098484718 0.5838784740291267 0 +-0.3774035668997982 0.56797066999209 0 +0.8562752679757717 0.3097481903629573 0 +-0.4557789974043656 -0.7431735260295338 0 +-0.4209719151177824 -0.8132782402637556 0 +-0.8751992713869907 -0.1008017497537349 0 +-0.8857226997807366 0.02999111817056395 0 +-0.7590491353497355 0.4406244770168601 0 +-0.791586086958344 0.09945532688856525 0 +0.4337503561672957 -0.7733007995086859 0 +0.690108180835496 -0.5609378978081738 0 +-0.7807171015479817 -0.4226616273600152 0 +-0.8766106165074704 -0.3056267462347799 0 +-0.5008014433296708 0.7269103600495153 0 +-0.4560818482234766 0.7994347002064576 0 +-0.3740194906333405 0.7732729413578556 0 +-0.1222511152728764 -0.8717183489764394 0 +0.9065022663437395 0.2522964138988663 0 +0.7630416073361792 -0.5181270503914466 0 +-0.3176700723976485 0.8297925006889542 0 +-0.2593031271325549 0.8896099725122354 0 +0.676915742067543 0.6223144879290337 0 +-0.68806459162302 -0.6360403687381796 0 +-0.7337914041536792 0.5724871772333908 0 +0.1101332462226518 -0.9305110966179267 0 +0.9330111819583282 0.1102505637103517 0 +0.5761266277257266 0.7385362993445825 0 +-0.5800947401331513 -0.7358466727989383 0 +0.248412640187938 -0.8878696047462332 0 +-0.2506998656961784 -0.9014056258200732 0 +0.8857712903621845 -0.2519961466192626 0 +0.5156180500524099 -0.7825912113462883 0 +0.9336126735260555 -0.09966181757432398 0 +0.04144716445137898 0.9381788894058519 0 +-0.1152218038618207 0.9372556142847692 0 +-0.8389514276272044 -0.1609622588059613 0 +-0.7095551930425135 0.4983839957254411 0 +0.7863894132866864 0.5265111215784929 0 +0.573902815104452 -0.5975726329253152 0 +-0.5513547726915167 0.669943413223456 0 +0.8555216260171135 -0.3944008703260372 0 +0.3245460811838435 -0.8857752720777493 0 +-0.2715530827423416 -0.8333831808528381 0 +0.8598241262673552 0.02259323489297163 0 +-0.9275666945403975 0.07294061119043314 0 +0.4631583886188241 0.8202659822638074 0 +-0.1967245053036725 -0.8523394114008389 0 +-0.6237586769749038 -0.6173908120224253 0 +-0.9391010515766621 -0.1111499343774911 0 +-0.3543964116638766 -0.8180390502763041 0 +0.597452793400188 0.6275506282918928 0 +-0.9108771185623138 0.259831626304533 0 +-0.5202157216983698 -0.7085383229879157 0 +-0.8254133743465917 -0.4591008553070199 0 +-0.8274017841490678 0.4556418125637993 0 +0.1676356747303709 -0.8574179890082834 0 +-0.5280699208077492 0.7867089745969889 0 +0.4626397746169288 -0.8286285533175972 0 +-0.9095868908081979 -0.2502579785736684 0 +0.5068862364459087 0.7115455276716114 0 +-0.105594754528827 -0.9396189004290079 0 +0.1637478285265344 0.8504535572789078 0 +0.2582205006795488 0.9155808353811785 0 +0.8513258167206499 -0.1874215809804928 0 +0.1118562598015134 0.945068770346427 0 +-0.3317948030005706 0.8900289912494354 0 +0.01343262898702855 0.8824404433752943 0 +-0.3879586989255454 -0.8656956926987225 0 +-0.1972982174631845 0.8550862030070706 0 +0.6805175931690028 0.5444269954909876 0 +0.6476605775476842 -0.70055190395823 0 +0.3149274546601034 0.8165819878834736 0 +0.6346180162227131 -0.6011285836236939 0 +-0.8270895021051141 0.04834245830966992 0 +0.9539770628486124 -0.03697069791150089 0 +-0.03755743158485551 -0.9559004687557859 0 +0.8243179538170899 -0.2944558788948262 0 +-0.9376721739981337 0.1271072706033549 0 +-0.9463792239571092 0.02111813918256253 0 +-0.8704203386385746 0.4012692708970406 0 +-0.7715646492774006 -0.3656789476977823 0 +0.5671355722288131 0.6828523091550884 0 +-0.7139113890002203 -0.5229602917555167 0 +0.8529708076007867 0.2306945719838603 0 +0.7448793710114532 -0.5911365238447537 0 +-0.4800779513212972 -0.8160765046231999 0 +$EndNodes +$Elements +5 1296 1 1296 +1 1 1 20 +1 1 5 +2 5 6 +3 6 7 +4 7 8 +5 8 9 +6 9 10 +7 10 11 +8 11 12 +9 12 13 +10 13 14 +11 14 15 +12 15 16 +13 16 17 +14 17 18 +15 18 19 +16 19 20 +17 20 21 +18 21 22 +19 22 23 +20 23 2 +1 2 1 20 +21 2 24 +22 24 25 +23 25 26 +24 26 27 +25 27 28 +26 28 29 +27 29 30 +28 30 31 +29 31 32 +30 32 33 +31 33 34 +32 34 35 +33 35 36 +34 36 37 +35 37 38 +36 38 39 +37 39 40 +38 40 41 +39 41 42 +40 42 3 +1 3 1 20 +41 3 43 +42 43 44 +43 44 45 +44 45 46 +45 46 47 +46 47 48 +47 48 49 +48 49 50 +49 50 51 +50 51 52 +51 52 53 +52 53 54 +53 54 55 +54 55 56 +55 56 57 +56 57 58 +57 58 59 +58 59 60 +59 60 61 +60 61 4 +1 4 1 20 +61 4 62 +62 62 63 +63 63 64 +64 64 65 +65 65 66 +66 66 67 +67 67 68 +68 68 69 +69 69 70 +70 70 71 +71 71 72 +72 72 73 +73 73 74 +74 74 75 +75 75 76 +76 76 77 +77 77 78 +78 78 79 +79 79 80 +80 80 1 +2 1 2 1216 +81 423 422 533 +82 8 9 566 +83 49 50 423 +84 85 479 570 +85 423 533 575 +86 373 87 374 +87 376 86 553 +88 479 85 569 +89 592 515 619 +90 50 99 423 +91 29 30 578 +92 73 74 582 +93 19 93 481 +94 372 88 589 +95 516 515 592 +96 553 86 598 +97 88 372 607 +98 47 90 576 +99 18 19 481 +100 102 365 599 +101 87 373 375 +102 365 102 622 +103 372 104 527 +104 516 91 530 +105 435 111 596 +106 504 92 518 +107 374 588 619 +108 414 108 415 +109 100 378 379 +110 111 435 627 +111 89 29 578 +112 9 92 566 +113 101 26 584 +114 104 372 589 +115 453 84 591 +116 588 103 619 +117 90 377 576 +118 100 379 595 +119 374 87 588 +120 103 592 619 +121 415 108 416 +122 369 98 633 +123 98 369 601 +124 8 566 581 +125 480 106 609 +126 518 92 519 +127 46 47 576 +128 106 480 623 +129 73 582 648 +130 569 102 599 +131 406 585 614 +132 74 112 582 +133 585 406 633 +134 92 504 566 +135 91 516 605 +136 379 573 595 +137 102 569 612 +138 115 587 600 +139 96 467 604 +140 530 91 573 +141 606 436 613 +142 467 96 640 +143 595 573 621 +144 109 405 625 +145 113 114 603 +146 26 27 584 +147 108 414 615 +148 405 109 628 +149 381 321 557 +150 382 112 467 +151 421 99 586 +152 578 30 620 +153 453 591 616 +154 49 423 617 +155 101 584 632 +156 95 606 613 +157 94 412 413 +158 378 100 380 +159 587 105 600 +160 375 373 381 +161 92 10 519 +162 10 11 519 +163 339 381 557 +164 373 321 381 +165 84 453 611 +166 481 93 635 +167 412 94 414 +168 590 623 645 +169 97 113 603 +170 369 518 519 +171 503 88 607 +172 98 585 633 +173 585 83 614 +174 415 416 572 +175 370 594 627 +176 594 370 640 +177 381 339 580 +178 516 592 605 +179 503 435 596 +180 580 339 610 +181 89 578 579 +182 331 412 414 +183 331 414 415 +184 384 119 564 +185 467 112 604 +186 340 384 564 +187 116 119 384 +188 331 415 417 +189 18 481 609 +190 117 383 386 +191 118 117 386 +192 252 270 417 +193 318 564 565 +194 270 331 417 +195 262 252 417 +196 389 387 560 +197 293 280 565 +198 318 340 564 +199 293 389 560 +200 123 387 389 +201 385 123 389 +202 120 123 385 +203 389 293 565 +204 121 118 386 +205 564 385 565 +206 280 318 565 +207 125 394 543 +208 121 386 563 +209 394 347 543 +210 394 121 563 +211 125 122 394 +212 291 347 563 +213 392 391 393 +214 390 391 392 +215 347 394 563 +216 128 131 390 +217 561 256 562 +218 393 391 558 +219 425 393 558 +220 390 131 391 +221 387 127 388 +222 127 128 390 +223 124 127 387 +224 388 127 390 +225 199 486 505 +226 486 359 505 +227 189 214 427 +228 230 200 399 +229 395 230 399 +230 200 199 399 +231 486 193 512 +232 214 426 427 +233 399 199 505 +234 359 486 512 +235 391 131 541 +236 337 296 404 +237 230 395 400 +238 297 401 402 +239 297 402 404 +240 193 427 512 +241 193 189 427 +242 337 404 405 +243 400 295 482 +244 296 297 404 +245 392 561 562 +246 31 32 97 +247 12 13 98 +248 50 51 99 +249 426 214 428 +250 227 398 550 +251 297 295 400 +252 230 400 482 +253 429 293 560 +254 407 428 439 +255 364 368 401 +256 462 179 469 +257 451 449 520 +258 69 70 100 +259 398 314 550 +260 408 273 452 +261 217 203 539 +262 297 400 401 +263 457 456 508 +264 409 408 452 +265 461 462 469 +266 272 273 410 +267 238 451 520 +268 409 452 453 +269 428 214 439 +270 388 390 392 +271 35 36 105 +272 16 17 106 +273 2 24 86 +274 195 457 508 +275 3 43 85 +276 273 408 410 +277 54 55 107 +278 395 364 401 +279 455 237 456 +280 4 62 87 +281 407 198 506 +282 1 5 88 +283 401 368 402 +284 430 411 493 +285 461 469 487 +286 398 227 551 +287 456 237 508 +288 424 217 539 +289 28 29 89 +290 272 410 411 +291 47 48 90 +292 396 125 543 +293 412 342 413 +294 179 462 538 +295 9 10 92 +296 488 179 538 +297 449 205 520 +298 25 26 101 +299 44 45 102 +300 292 342 412 +301 162 170 463 +302 463 461 487 +303 400 395 401 +304 237 455 458 +305 330 353 354 +306 324 407 506 +307 6 7 104 +308 170 208 463 +309 198 407 439 +310 554 397 556 +311 407 326 428 +312 366 466 467 +313 366 346 466 +314 431 432 433 +315 63 64 103 +316 411 410 493 +317 205 449 450 +318 411 430 437 +319 443 554 556 +320 532 309 533 +321 466 346 498 +322 212 206 450 +323 66 67 91 +324 356 354 357 +325 418 420 452 +326 422 532 533 +327 452 420 453 +328 328 329 330 +329 307 309 528 +330 418 419 420 +331 437 164 485 +332 431 253 432 +333 249 293 429 +334 184 186 458 +335 292 341 342 +336 426 355 427 +337 402 368 403 +338 411 437 485 +339 438 458 492 +340 433 432 537 +341 330 354 355 +342 330 355 426 +343 206 205 450 +344 272 411 485 +345 529 217 534 +346 296 295 297 +347 289 285 332 +348 202 230 482 +349 204 202 482 +350 197 229 233 +351 330 329 353 +352 278 434 435 +353 126 125 396 +354 326 327 328 +355 457 195 529 +356 448 544 545 +357 331 292 412 +358 328 330 426 +359 321 314 398 +360 254 253 431 +361 197 233 272 +362 355 354 356 +363 188 187 189 +364 348 498 522 +365 313 431 433 +366 217 424 534 +367 328 327 329 +368 528 309 532 +369 289 332 382 +370 358 357 476 +371 424 350 534 +372 265 349 460 +373 367 480 481 +374 270 291 292 +375 324 326 407 +376 235 265 460 +377 246 228 247 +378 448 451 544 +379 252 269 270 +380 282 247 284 +381 240 251 252 +382 261 254 431 +383 318 319 340 +384 226 191 228 +385 361 358 476 +386 458 455 492 +387 186 237 458 +388 287 313 314 +389 272 233 273 +390 313 261 431 +391 197 272 485 +392 180 177 488 +393 287 261 313 +394 344 408 409 +395 235 263 264 +396 265 264 266 +397 250 178 438 +398 160 164 437 +399 250 438 492 +400 159 155 160 +401 356 357 358 +402 354 353 406 +403 163 276 474 +404 168 172 197 +405 324 325 326 +406 186 185 187 +407 235 264 265 +408 498 346 522 +409 284 343 344 +410 186 187 188 +411 241 212 450 +412 457 529 534 +413 284 345 436 +414 228 191 430 +415 270 292 331 +416 215 216 440 +417 290 346 366 +418 266 285 289 +419 326 328 428 +420 438 184 458 +421 342 341 383 +422 221 239 240 +423 171 163 474 +424 191 166 430 +425 216 225 440 +426 131 132 541 +427 162 463 487 +428 284 344 345 +429 249 280 293 +430 445 444 446 +431 268 274 290 +432 320 314 321 +433 360 358 361 +434 177 179 488 +435 367 362 480 +436 294 258 464 +437 454 455 456 +438 277 239 442 +439 257 222 258 +440 223 204 482 +441 362 361 480 +442 203 215 539 +443 460 349 461 +444 246 247 282 +445 225 256 440 +446 166 160 437 +447 159 160 166 +448 236 249 429 +449 258 222 464 +450 190 189 193 +451 443 444 445 +452 284 247 343 +453 178 181 438 +454 252 251 269 +455 304 323 324 +456 154 155 159 +457 326 325 327 +458 240 239 251 +459 208 234 235 +460 181 183 184 +461 430 166 437 +462 226 228 246 +463 352 454 456 +464 225 236 256 +465 432 474 475 +466 239 220 442 +467 282 284 436 +468 474 276 475 +469 273 418 452 +470 259 223 295 +471 328 426 428 +472 318 317 319 +473 290 370 434 +474 290 366 370 +475 290 274 346 +476 270 269 291 +477 154 151 155 +478 445 446 447 +479 221 220 239 +480 352 275 454 +481 135 136 490 +482 259 295 296 +483 443 241 444 +484 461 349 462 +485 442 241 443 +486 233 267 418 +487 266 264 285 +488 136 139 490 +489 352 456 457 +490 228 430 493 +491 265 266 348 +492 208 460 463 +493 465 135 490 +494 257 258 268 +495 344 343 408 +496 222 218 464 +497 208 235 460 +498 364 362 367 +499 256 236 429 +500 446 449 451 +501 338 282 436 +502 292 291 341 +503 280 317 318 +504 220 241 442 +505 160 155 161 +506 273 233 418 +507 218 180 488 +508 210 307 494 +509 268 258 274 +510 447 446 448 +511 268 290 434 +512 132 134 135 +513 184 185 186 +514 188 189 190 +515 466 382 467 +516 138 137 472 +517 304 324 506 +518 128 130 131 +519 235 234 263 +520 132 133 134 +521 446 444 449 +522 135 134 136 +523 33 34 82 +524 220 212 241 +525 132 130 133 +526 184 183 185 +527 219 218 222 +528 298 147 468 +529 129 397 459 +530 249 279 280 +531 213 220 221 +532 343 247 493 +533 275 271 454 +534 357 354 406 +535 181 184 438 +536 454 301 455 +537 136 134 137 +538 253 171 474 +539 141 138 472 +540 448 446 451 +541 136 137 138 +542 432 253 474 +543 245 222 257 +544 194 193 486 +545 285 264 497 +546 147 144 468 +547 271 301 454 +548 153 156 157 +549 274 258 294 +550 199 194 486 +551 176 178 250 +552 131 130 132 +553 144 141 468 +554 179 175 469 +555 300 465 490 +556 146 144 147 +557 242 275 300 +558 324 323 325 +559 295 223 482 +560 170 207 208 +561 168 197 485 +562 198 491 506 +563 150 151 154 +564 468 141 472 +565 151 147 298 +566 161 155 496 +567 309 308 521 +568 133 130 459 +569 146 147 150 +570 142 141 144 +571 160 161 164 +572 271 299 301 +573 236 248 249 +574 300 275 350 +575 140 138 141 +576 149 153 441 +577 350 275 352 +578 300 351 465 +579 238 161 496 +580 153 174 441 +581 150 147 151 +582 420 419 421 +583 242 271 275 +584 14 15 83 +585 186 188 237 +586 197 172 229 +587 210 232 307 +588 175 157 469 +589 145 149 483 +590 441 507 513 +591 307 308 309 +592 189 187 214 +593 360 361 362 +594 142 145 483 +595 181 182 183 +596 278 268 434 +597 470 262 478 +598 178 180 181 +599 153 157 174 +600 174 175 176 +601 196 194 199 +602 257 268 278 +603 167 166 191 +604 202 200 230 +605 348 266 498 +606 251 239 277 +607 269 251 495 +608 139 242 490 +609 418 267 419 +610 174 157 175 +611 507 473 513 +612 548 227 550 +613 242 143 271 +614 537 475 542 +615 286 261 287 +616 347 269 495 +617 283 282 338 +618 300 350 351 +619 410 343 493 +620 233 229 267 +621 140 141 142 +622 289 466 498 +623 271 143 299 +624 240 252 262 +625 301 299 473 +626 149 152 153 +627 145 144 146 +628 153 152 156 +629 289 382 466 +630 201 203 217 +631 247 228 493 +632 177 175 179 +633 364 367 368 +634 265 348 349 +635 164 168 485 +636 444 241 450 +637 433 548 550 +638 176 177 178 +639 260 254 261 +640 225 231 236 +641 349 348 535 +642 351 350 424 +643 280 279 317 +644 301 473 492 +645 52 53 84 +646 216 204 223 +647 221 240 470 +648 491 302 506 +649 460 461 463 +650 451 238 544 +651 130 129 459 +652 475 511 542 +653 163 158 276 +654 266 289 498 +655 149 441 483 +656 455 301 492 +657 185 198 439 +658 545 472 547 +659 351 424 425 +660 128 129 130 +661 473 250 492 +662 150 154 306 +663 213 212 220 +664 208 207 234 +665 139 138 140 +666 441 174 507 +667 190 193 194 +668 151 298 496 +669 157 156 487 +670 469 157 487 +671 287 314 320 +672 142 144 145 +673 227 226 551 +674 129 396 397 +675 139 143 242 +676 209 171 253 +677 332 285 371 +678 320 321 373 +679 146 150 477 +680 187 185 439 +681 213 221 244 +682 176 250 507 +683 136 138 139 +684 242 300 490 +685 475 276 511 +686 142 483 499 +687 211 206 212 +688 249 248 279 +689 464 218 488 +690 129 126 396 +691 182 218 219 +692 250 473 507 +693 276 158 484 +694 359 358 360 +695 145 148 149 +696 182 180 218 +697 176 175 177 +698 278 435 502 +699 155 151 496 +700 192 191 226 +701 159 166 167 +702 145 146 148 +703 140 142 499 +704 156 162 487 +705 190 195 508 +706 291 269 347 +707 449 444 450 +708 149 148 152 +709 139 140 143 +710 274 294 522 +711 371 285 497 +712 471 470 478 +713 236 231 248 +714 164 165 168 +715 264 263 497 +716 312 296 337 +717 148 477 484 +718 244 221 470 +719 165 161 238 +720 128 126 129 +721 251 277 495 +722 263 305 315 +723 196 199 200 +724 148 146 477 +725 188 190 508 +726 164 161 165 +727 150 306 477 +728 178 177 180 +729 203 204 215 +730 211 212 213 +731 156 158 162 +732 185 183 198 +733 307 232 308 +734 244 470 471 +735 259 296 312 +736 201 202 203 +737 336 303 489 +738 408 343 410 +739 209 253 254 +740 183 182 491 +741 502 435 503 +742 190 194 195 +743 299 143 499 +744 224 223 259 +745 156 152 158 +746 152 148 484 +747 246 282 283 +748 260 261 286 +749 304 322 323 +750 162 163 170 +751 170 171 207 +752 237 188 508 +753 181 180 182 +754 308 311 521 +755 263 315 497 +756 263 234 305 +757 311 365 377 +758 299 499 513 +759 195 194 196 +760 133 459 552 +761 196 200 201 +762 219 222 245 +763 473 299 513 +764 174 176 507 +765 302 303 304 +766 216 224 225 +767 168 169 172 +768 501 244 517 +769 216 223 224 +770 240 262 470 +771 356 358 359 +772 225 224 231 +773 173 169 205 +774 201 217 529 +775 127 126 128 +776 219 302 491 +777 201 200 202 +778 499 483 513 +779 143 140 499 +780 363 362 364 +781 173 205 206 +782 327 325 504 +783 245 257 489 +784 315 316 378 +785 182 219 491 +786 210 206 211 +787 427 355 512 +788 203 202 204 +789 215 204 216 +790 477 306 511 +791 162 158 163 +792 170 163 171 +793 173 206 210 +794 168 165 169 +795 302 304 506 +796 283 338 339 +797 383 341 386 +798 210 211 232 +799 317 279 335 +800 211 213 501 +801 306 154 523 +802 198 183 491 +803 483 441 513 +804 158 152 484 +805 248 231 281 +806 484 477 511 +807 350 352 534 +808 207 171 209 +809 308 310 311 +810 154 159 523 +811 287 320 525 +812 363 395 399 +813 255 254 260 +814 310 501 517 +815 255 260 288 +816 219 245 302 +817 209 254 255 +818 308 232 310 +819 322 303 336 +820 311 310 334 +821 234 207 500 +822 552 459 554 +823 360 362 363 +824 346 274 522 +825 464 488 538 +826 159 167 523 +827 311 377 521 +828 257 278 489 +829 355 356 512 +830 243 537 542 +831 514 515 516 +832 167 191 192 +833 303 245 489 +834 172 169 173 +835 232 211 501 +836 213 244 501 +837 209 255 500 +838 315 305 316 +839 302 245 303 +840 310 232 501 +841 304 303 322 +842 279 248 509 +843 359 360 505 +844 165 238 520 +845 363 364 395 +846 244 471 517 +847 167 192 243 +848 459 397 554 +849 229 172 494 +850 207 209 500 +851 325 323 333 +852 353 329 369 +853 276 484 511 +854 305 234 500 +855 195 196 529 +856 286 287 525 +857 462 349 535 +858 322 336 372 +859 471 478 479 +860 167 243 523 +861 214 187 439 +862 224 259 524 +863 173 210 494 +864 336 489 502 +865 356 359 512 +866 348 522 535 +867 514 286 515 +868 334 310 517 +869 320 374 525 +870 329 327 518 +871 205 169 520 +872 172 173 494 +873 314 313 550 +874 320 373 374 +875 494 307 528 +876 288 260 514 +877 259 312 524 +878 260 286 514 +879 231 224 524 +880 360 363 505 +881 192 226 227 +882 248 281 509 +883 132 135 541 +884 335 279 509 +885 311 334 365 +886 255 288 531 +887 363 399 505 +888 169 165 520 +889 312 337 376 +890 335 509 510 +891 468 472 545 +892 281 231 524 +893 196 201 529 +894 515 286 525 +895 489 278 502 +896 40 41 108 +897 323 322 527 +898 535 294 538 +899 21 22 109 +900 511 306 542 +901 462 535 538 +902 352 457 534 +903 522 294 535 +904 325 333 504 +905 327 504 518 +906 378 316 379 +907 238 496 544 +908 421 419 422 +909 333 323 527 +910 59 60 110 +911 517 471 526 +912 393 425 559 +913 267 229 528 +914 38 39 94 +915 74 75 112 +916 448 545 547 +917 19 20 93 +918 432 475 537 +919 334 517 526 +920 288 514 530 +921 514 516 530 +922 294 464 538 +923 267 528 532 +924 135 465 541 +925 309 521 533 +926 134 133 549 +927 229 494 528 +928 57 58 95 +929 316 305 531 +930 442 443 556 +931 305 500 531 +932 322 372 527 +933 471 479 526 +934 425 424 539 +935 369 329 518 +936 124 126 127 +937 347 495 543 +938 500 255 531 +939 78 79 111 +940 288 530 546 +941 419 267 532 +942 544 298 545 +943 433 537 548 +944 531 288 546 +945 226 246 551 +946 306 523 542 +947 509 281 536 +948 510 509 536 +949 496 298 544 +950 422 419 532 +951 524 312 540 +952 536 281 540 +953 298 468 545 +954 312 376 540 +955 76 77 96 +956 243 192 548 +957 316 531 546 +958 447 448 547 +959 397 396 555 +960 523 243 542 +961 472 137 547 +962 447 547 549 +963 281 524 540 +964 397 555 556 +965 547 137 549 +966 537 243 548 +967 549 133 552 +968 543 495 555 +969 137 134 549 +970 396 543 555 +971 313 433 550 +972 379 316 546 +973 246 283 551 +974 192 227 548 +975 445 447 552 +976 321 398 557 +977 71 72 81 +978 445 552 554 +979 443 445 554 +980 277 442 556 +981 536 540 553 +982 495 277 555 +983 447 549 552 +984 555 277 556 +985 540 376 553 +986 351 425 558 +987 541 465 558 +988 391 541 558 +989 551 283 557 +990 124 125 126 +991 398 551 557 +992 283 339 557 +993 256 429 562 +994 559 440 561 +995 465 351 558 +996 425 539 559 +997 215 440 559 +998 393 559 561 +999 392 393 561 +1000 123 124 387 +1001 124 122 125 +1002 539 215 559 +1003 429 560 562 +1004 387 388 560 +1005 440 256 561 +1006 388 392 562 +1007 560 388 562 +1008 123 122 124 +1009 341 291 563 +1010 122 121 394 +1011 120 122 123 +1012 386 341 563 +1013 120 121 122 +1014 120 118 121 +1015 119 120 385 +1016 119 385 564 +1017 385 389 565 +1018 119 118 120 +1019 116 118 119 +1020 116 117 118 +1021 116 115 117 +1022 114 116 384 +1023 504 333 566 +1024 344 409 567 +1025 345 567 568 +1026 345 344 567 +1027 114 115 116 +1028 342 383 571 +1029 478 262 572 +1030 579 335 583 +1031 417 415 572 +1032 526 479 569 +1033 413 342 571 +1034 262 417 572 +1035 114 82 115 +1036 546 530 573 +1037 527 104 647 +1038 403 109 625 +1039 574 371 636 +1040 81 574 636 +1041 88 503 638 +1042 317 335 579 +1043 582 574 648 +1044 423 99 646 +1045 422 423 646 +1046 101 553 598 +1047 578 319 579 +1048 89 579 583 +1049 112 382 582 +1050 379 546 573 +1051 421 586 611 +1052 382 332 582 +1053 606 593 610 +1054 553 101 632 +1055 577 319 578 +1056 590 106 623 +1057 332 371 574 +1058 115 82 587 +1059 571 383 600 +1060 377 365 576 +1061 378 380 602 +1062 340 319 577 +1063 96 594 640 +1064 117 115 600 +1065 365 334 599 +1066 105 571 600 +1067 570 479 637 +1068 526 569 599 +1069 319 317 579 +1070 110 580 610 +1071 375 381 580 +1072 109 403 626 +1073 114 384 603 +1074 623 476 645 +1075 315 378 602 +1076 371 497 602 +1077 481 480 609 +1078 377 90 644 +1079 113 82 114 +1080 593 110 610 +1081 104 7 581 +1082 7 8 581 +1083 594 111 627 +1084 56 568 649 +1085 372 336 607 +1086 502 503 607 +1087 568 567 649 +1088 567 107 649 +1089 87 375 639 +1090 413 571 618 +1091 93 403 635 +1092 436 345 613 +1093 476 357 614 +1094 576 365 622 +1095 335 510 583 +1096 525 374 619 +1097 32 33 113 +1098 13 14 585 +1099 51 52 586 +1100 33 82 113 +1101 14 83 585 +1102 52 84 586 +1103 97 32 113 +1104 98 13 585 +1105 99 51 586 +1106 6 104 589 +1107 34 35 587 +1108 82 34 587 +1109 63 103 588 +1110 15 16 590 +1111 35 105 587 +1112 103 64 592 +1113 53 54 591 +1114 83 15 590 +1115 16 106 590 +1116 84 53 591 +1117 54 107 591 +1118 62 63 588 +1119 5 6 589 +1120 64 65 592 +1121 87 62 588 +1122 88 5 589 +1123 58 59 593 +1124 77 78 594 +1125 95 58 593 +1126 59 110 593 +1127 96 77 594 +1128 78 111 594 +1129 453 420 611 +1130 375 580 624 +1131 369 519 601 +1132 69 100 595 +1133 79 80 596 +1134 23 2 597 +1135 480 361 623 +1136 405 404 625 +1137 402 403 625 +1138 479 478 637 +1139 68 69 595 +1140 111 79 596 +1141 435 434 627 +1142 2 86 597 +1143 337 405 630 +1144 577 97 603 +1145 340 577 603 +1146 536 553 632 +1147 25 101 598 +1148 86 24 598 +1149 533 521 644 +1150 353 369 633 +1151 566 333 647 +1152 366 467 640 +1153 367 481 635 +1154 403 368 635 +1155 332 574 582 +1156 416 570 637 +1157 86 376 630 +1158 94 413 643 +1159 519 11 601 +1160 11 12 601 +1161 112 75 604 +1162 76 96 604 +1163 567 409 616 +1164 24 25 598 +1165 375 624 639 +1166 66 91 605 +1167 107 567 616 +1168 334 526 599 +1169 383 117 600 +1170 575 533 644 +1171 403 93 626 +1172 17 18 609 +1173 570 416 608 +1174 414 94 615 +1175 583 510 584 +1176 12 98 601 +1177 44 102 612 +1178 85 43 612 +1179 581 566 647 +1180 568 56 631 +1181 94 39 615 +1182 40 108 615 +1183 90 48 617 +1184 105 36 618 +1185 569 85 612 +1186 380 100 634 +1187 81 380 634 +1188 571 105 618 +1189 31 97 620 +1190 91 67 621 +1191 45 46 622 +1192 110 60 624 +1193 93 20 626 +1194 21 109 626 +1195 338 436 606 +1196 109 22 628 +1197 405 597 630 +1198 497 315 602 +1199 413 618 643 +1200 577 578 620 +1201 55 56 649 +1202 81 72 648 +1203 100 70 634 +1204 71 81 634 +1205 38 94 643 +1206 42 3 642 +1207 1 88 638 +1208 41 42 641 +1209 4 87 639 +1210 384 340 603 +1211 584 510 632 +1212 75 76 604 +1213 575 90 617 +1214 65 66 605 +1215 573 91 621 +1216 572 416 637 +1217 336 502 607 +1218 46 576 622 +1219 97 577 620 +1220 99 421 646 +1221 574 81 648 +1222 592 65 605 +1223 57 95 631 +1224 85 570 642 +1225 416 108 641 +1226 580 110 624 +1227 106 17 609 +1228 28 89 629 +1229 409 453 616 +1230 339 338 610 +1231 90 575 644 +1232 420 421 611 +1233 43 44 612 +1234 423 575 617 +1235 345 568 613 +1236 23 597 628 +1237 104 581 647 +1238 357 406 614 +1239 338 606 610 +1240 39 40 615 +1241 68 595 621 +1242 597 405 628 +1243 48 49 617 +1244 596 80 638 +1245 36 37 618 +1246 503 596 638 +1247 515 525 619 +1248 95 593 606 +1249 421 422 646 +1250 30 31 620 +1251 67 68 621 +1252 102 45 622 +1253 361 476 623 +1254 89 583 629 +1255 27 28 629 +1256 608 42 642 +1257 42 608 641 +1258 60 61 624 +1259 404 402 625 +1260 586 84 611 +1261 56 57 631 +1262 20 21 626 +1263 434 370 627 +1264 22 23 628 +1265 380 81 636 +1266 476 614 645 +1267 376 337 630 +1268 371 602 636 +1269 510 536 632 +1270 406 353 633 +1271 591 107 616 +1272 70 71 634 +1273 368 367 635 +1274 478 572 637 +1275 80 1 638 +1276 61 4 639 +1277 370 366 640 +1278 108 41 641 +1279 3 85 642 +1280 37 38 643 +1281 521 377 644 +1282 333 527 647 +1283 72 73 648 +1284 107 55 649 +1285 95 613 631 +1286 584 27 629 +1287 618 37 643 +1288 597 86 630 +1289 624 61 639 +1290 583 584 629 +1291 570 608 642 +1292 608 416 641 +1293 83 590 645 +1294 602 380 636 +1295 613 568 631 +1296 614 83 645 +$EndElements diff --git a/models/lshape.msh b/models/lshape.msh new file mode 100644 index 0000000..ab453ce --- /dev/null +++ b/models/lshape.msh @@ -0,0 +1,1680 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +2 +1 1 "boundary" +2 2 "domain" +$EndPhysicalNames +$Entities +6 6 1 0 +1 0 0 0 0 +2 1 0 0 0 +3 1 0.5 0 0 +4 0.5 0.5 0 0 +5 0.5 1 0 0 +6 0 1 0 0 +1 0 0 0 1 0 0 1 1 2 1 -2 +2 1 0 0 1 0.5 0 1 1 2 2 -3 +3 0.5 0.5 0 1 0.5 0 1 1 2 3 -4 +4 0.5 0.5 0 0.5 1 0 1 1 2 4 -5 +5 0 1 0 0.5 1 0 1 1 2 5 -6 +6 0 0 0 0 1 0 1 1 2 6 -1 +1 0 0 0 1 1 0 1 2 6 1 2 3 4 5 6 +$EndEntities +$Nodes +13 408 1 408 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +1 0 0 +0 3 0 1 +3 +1 0.5 0 +0 4 0 1 +4 +0.5 0.5 0 +0 5 0 1 +5 +0.5 1 0 +0 6 0 1 +6 +0 1 0 +1 1 0 19 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +0.04999999999989967 0 0 +0.09999999999981468 0 0 +0.1499999999997036 0 0 +0.1999999999995579 0 0 +0.2499999999994122 0 0 +0.2999999999992664 0 0 +0.3499999999991207 0 0 +0.399999999998975 0 0 +0.4499999999988292 0 0 +0.4999999999986943 0 0 +0.5499999999988151 0 0 +0.5999999999989468 0 0 +0.6499999999990784 0 0 +0.69999999999921 0 0 +0.7499999999993417 0 0 +0.7999999999994734 0 0 +0.8499999999996051 0 0 +0.8999999999997368 0 0 +0.9499999999998684 0 0 +1 2 0 9 +26 +27 +28 +29 +30 +31 +32 +33 +34 +1 0.04999999999990734 0 +1 0.09999999999977895 0 +1 0.1499999999996332 0 +1 0.1999999999994875 0 +1 0.2499999999993471 0 +1 0.2999999999994734 0 +1 0.349999999999605 0 +1 0.3999999999997367 0 +1 0.4499999999998684 0 +1 3 0 9 +35 +36 +37 +38 +39 +40 +41 +42 +43 +0.9500000000000001 0.5 0 +0.9000000000000002 0.5 0 +0.8500000000000002 0.5 0 +0.8000000000000002 0.5 0 +0.7500000000000001 0.5 0 +0.7000000000000001 0.5 0 +0.6500000000000001 0.5 0 +0.6000000000000001 0.5 0 +0.55 0.5 0 +1 4 0 9 +44 +45 +46 +47 +48 +49 +50 +51 +52 +0.5 0.5499999999999999 0 +0.5 0.5999999999999998 0 +0.5 0.6499999999999998 0 +0.5 0.6999999999999998 0 +0.5 0.7499999999999999 0 +0.5 0.7999999999999999 0 +0.5 0.8499999999999999 0 +0.5 0.8999999999999999 0 +0.5 0.95 0 +1 5 0 9 +53 +54 +55 +56 +57 +58 +59 +60 +61 +0.4499999999997918 1 0 +0.3999999999999999 1 0 +0.3500000000003467 1 0 +0.3000000000006934 1 0 +0.2500000000010293 1 0 +0.2000000000008322 1 0 +0.1500000000006241 1 0 +0.1000000000004161 1 0 +0.05000000000020799 1 0 +1 6 0 19 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +0 0.9499999999997918 0 +0 0.8999999999995836 0 +0 0.8499999999996529 0 +0 0.7999999999999998 0 +0 0.7500000000003466 0 +0 0.7000000000006934 0 +0 0.6500000000010401 0 +0 0.6000000000013869 0 +0 0.5500000000017335 0 +0 0.5000000000020587 0 +0 0.4500000000018723 0 +0 0.4000000000016644 0 +0 0.3500000000014564 0 +0 0.3000000000012483 0 +0 0.2500000000010403 0 +0 0.2000000000008322 0 +0 0.1500000000006241 0 +0 0.100000000000416 0 +0 0.05000000000020799 0 +2 1 0 328 +81 +82 +83 +84 +85 +86 +87 +88 +89 +90 +91 +92 +93 +94 +95 +96 +97 +98 +99 +100 +101 +102 +103 +104 +105 +106 +107 +108 +109 +110 +111 +112 +113 +114 +115 +116 +117 +118 +119 +120 +121 +122 +123 +124 +125 +126 +127 +128 +129 +130 +131 +132 +133 +134 +135 +136 +137 +138 +139 +140 +141 +142 +143 +144 +145 +146 +147 +148 +149 +150 +151 +152 +153 +154 +155 +156 +157 +158 +159 +160 +161 +162 +163 +164 +165 +166 +167 +168 +169 +170 +171 +172 +173 +174 +175 +176 +177 +178 +179 +180 +181 +182 +183 +184 +185 +186 +187 +188 +189 +190 +191 +192 +193 +194 +195 +196 +197 +198 +199 +200 +201 +202 +203 +204 +205 +206 +207 +208 +209 +210 +211 +212 +213 +214 +215 +216 +217 +218 +219 +220 +221 +222 +223 +224 +225 +226 +227 +228 +229 +230 +231 +232 +233 +234 +235 +236 +237 +238 +239 +240 +241 +242 +243 +244 +245 +246 +247 +248 +249 +250 +251 +252 +253 +254 +255 +256 +257 +258 +259 +260 +261 +262 +263 +264 +265 +266 +267 +268 +269 +270 +271 +272 +273 +274 +275 +276 +277 +278 +279 +280 +281 +282 +283 +284 +285 +286 +287 +288 +289 +290 +291 +292 +293 +294 +295 +296 +297 +298 +299 +300 +301 +302 +303 +304 +305 +306 +307 +308 +309 +310 +311 +312 +313 +314 +315 +316 +317 +318 +319 +320 +321 +322 +323 +324 +325 +326 +327 +328 +329 +330 +331 +332 +333 +334 +335 +336 +337 +338 +339 +340 +341 +342 +343 +344 +345 +346 +347 +348 +349 +350 +351 +352 +353 +354 +355 +356 +357 +358 +359 +360 +361 +362 +363 +364 +365 +366 +367 +368 +369 +370 +371 +372 +373 +374 +375 +376 +377 +378 +379 +380 +381 +382 +383 +384 +385 +386 +387 +388 +389 +390 +391 +392 +393 +394 +395 +396 +397 +398 +399 +400 +401 +402 +403 +404 +405 +406 +407 +408 +0.3273509456663035 0.04703283710292606 0 +0.03913125464393551 0.3241179869740647 0 +0.04226795090231397 0.7253214333193417 0 +0.04226795090225829 0.5746785666825617 0 +0.4749999999987617 0.04330127018910039 0 +0.7246594085730371 0.4575893408509423 0 +0.5769337567297425 0.4578151847792432 0 +0.4566987298107837 0.5249999999999877 0 +0.4566987298107846 0.4750000000000113 0 +0.4566987298107066 0.6749999999999028 0 +0.6749999999989884 0.04330127018955592 0 +0.413397459621586 0.4999999999999999 0 +0.4133974596215988 0.4500000000000464 0 +0.2734201453150389 0.957610859339088 0 +0.956698729810902 0.2249999999994173 0 +0.3700961894324296 0.475000000000022 0 +0.3700961894324335 0.425000000000088 0 +0.4133974596215887 0.4000000000001159 0 +0.3700961894324262 0.375000000000168 0 +0.3267949192432816 0.4500000000000464 0 +0.3267949192432613 0.4999999999999659 0 +0.413397459621567 0.3500000000002013 0 +0.04330127018940215 0.1750000000007281 0 +0.3700961894324188 0.3250000000002572 0 +0.2834936490541388 0.4749999999999909 0 +0.2834936490541207 0.5249999999998722 0 +0.4133974596215535 0.3000000000002928 0 +0.4566987298107781 0.825 0 +0.8249999999995391 0.04330127018933595 0 +0.3700961894324066 0.2750000000003497 0 +0.3267949192432658 0.300000000000319 0 +0.1745399000490335 0.04442715103593312 0 +0.3267949192432517 0.250000000000393 0 +0.2834936490540974 0.275000000000373 0 +0.4566987298107036 0.3250000000002339 0 +0.4566987298106964 0.2750000000003251 0 +0.499999999999854 0.3000000000002635 0 +0.4999999999998446 0.2500000000003421 0 +0.5423891406608721 0.2734201453144625 0 +0.5432090792423692 0.2248403205967806 0 +0.5855230241042531 0.2481302229662452 0 +0.2401923788650244 0.4999999999999259 0 +0.2401923788650047 0.5499999999997811 0 +0.1750000000007831 0.9566987298106624 0 +0.9566987298106664 0.3249999999995392 0 +0.2834936490540848 0.5749999999997262 0 +0.2401923788649808 0.5999999999996251 0 +0.2834936490540436 0.6249999999995599 0 +0.2401923788649501 0.6499999999994512 0 +0.2834936490540146 0.6749999999993914 0 +0.2401923788649157 0.6999999999992831 0 +0.2834936490539851 0.7249999999992316 0 +0.2401923788648786 0.7499999999991258 0 +0.04330127018939268 0.4250000000017683 0 +0.8257265603118509 0.4628218759961584 0 +0.04258834754699556 0.826234818237667 0 +0.2401923788650402 0.4500000000000848 0 +0.2834936490540791 0.2250000000004266 0 +0.2401923788649036 0.2500000000004211 0 +0.3700961894324005 0.2250000000004342 0 +0.2834936490539628 0.7749999999990829 0 +0.5749999999987638 0.0433012701893111 0 +0.5881815682283751 0.202225884689311 0 +0.6299038105673903 0.2250000000003173 0 +0.628135280138798 0.2800495745688839 0 +0.6735736926147347 0.2496926105356483 0 +0.2401923788648558 0.7999999999989853 0 +0.2408903188277625 0.2012088674767907 0 +0.3793362827355576 0.955185193752226 0 +0.9602038715518374 0.1270105983836034 0 +0.6719519356083798 0.2031074184742508 0 +0.7161692731969996 0.2257951645223003 0 +0.7171749735602573 0.2739323106170902 0 +0.7598076211349941 0.2500000000002988 0 +0.7598076211349929 0.2000000000003314 0 +0.760292743529404 0.2991748098579431 0 +0.1968911086759257 0.5249999999998386 0 +0.196891108675821 0.6749999999993914 0 +0.1968911086757669 0.7749999999990328 0 +0.8031897450565939 0.2748624683098995 0 +0.197720354644937 0.2264362961511769 0 +0.1970293163372437 0.2752393826922421 0 +0.4992071242570729 0.2015273340245827 0 +0.413397459621584 0.549999999999961 0 +0.3267949192430616 0.7499999999991684 0 +0.3267949192430466 0.7999999999990255 0 +0.1544640034003508 0.2515140980455522 0 +0.7099023383753856 0.3265847429648525 0 +0.2838794906441409 0.1756682972381371 0 +0.5474303647033962 0.3265474465536659 0 +0.1963491463203765 0.8240612936635997 0 +0.2046079674708329 0.1763939268543614 0 +0.8035769059706621 0.3241923044496866 0 +0.326794919243123 0.6499999999995505 0 +0.283493649054134 0.3250000000003258 0 +0.1537585672490432 0.3002922467899907 0 +0.1551466452653616 0.7998435489430606 0 +0.1535898384866366 0.8499999999987833 0 +0.846501639576541 0.2998424621267869 0 +0.8464176380767888 0.2499871984955454 0 +0.8465034102983295 0.3498391277629294 0 +0.2373988665767464 0.8463154102742384 0 +0.04789563223317055 0.9206942515958663 0 +0.04218481522090629 0.07693375673009624 0 +0.9614259745565655 0.424769995604611 0 +0.3700961894321509 0.7749999999991137 0 +0.1538493062831467 0.7499739248231334 0 +0.37009618943212 0.8249999999989704 0 +0.3267949192430335 0.8499999999988717 0 +0.1969044250206779 0.3250230645864978 0 +0.1536201793377592 0.3500525518964682 0 +0.3700961894321292 0.7249999999993213 0 +0.8041534681084057 0.3705581194213567 0 +0.5000000000000011 0.4461324865405391 0 +0.4576613977472225 0.9247881315359677 0 +0.9247033841501502 0.04147487020751377 0 +0.1111753065519996 0.2740662392349797 0 +0.114674807006437 0.7758602837129869 0 +0.8897114317025984 0.2750000000002901 0 +0.413397459621346 0.7499999999994161 0 +0.1968983848751689 0.3750126027474381 0 +0.1535961079950131 0.4000108591076815 0 +0.1102885682972755 0.3750000000006132 0 +0.1171204662489031 0.2250967228804339 0 +0.3692588482059132 0.873549682451865 0 +0.3266900717439968 0.8984027906278738 0 +0.1117667103449756 0.7250506686422695 0 +0.8897114317026018 0.2250000000003171 0 +0.8464114076073088 0.1999978664162154 0 +0.8897114317026023 0.1750000000003673 0 +0.8464103691957262 0.1499996444030262 0 +0.8891449553730094 0.125946596873505 0 +0.1102885682973172 0.4250000000006584 0 +0.1535908834047255 0.4500018098516126 0 +0.1102885682973937 0.4750000000006693 0 +0.07572130213083754 0.2461560389724615 0 +0.8474065172805127 0.3959009913615256 0 +0.88971143170262 0.3750000000002827 0 +0.8857803986017818 0.4219143686855498 0 +0.8483806022701861 0.09270099726907755 0 +0.8022283194375899 0.1224479899637018 0 +0.04330127018916183 0.4750000000019656 0 +0.2247366908866135 0.9568507513988467 0 +0.2034335111761244 0.9156959280729696 0 +0.1529655533773328 0.9054135580583327 0 +0.2504716516569284 0.9180818635997836 0 +0.1250000000007554 0.9566987298104703 0 +0.1002204983512368 0.9169846130835413 0 +0.1538794401273502 0.7000040989105001 0 +0.1122303260836102 0.6750091279251574 0 +0.1539617317245088 0.6500022044723518 0 +0.1123213106394492 0.6250018887327804 0 +0.1539906110833617 0.6000006822007331 0 +0.1112484087764576 0.575000111539232 0 +0.07605655127800627 0.6500018361096801 0 +0.1538166073326923 0.5500001322900299 0 +0.04330127018915887 0.5250000000017628 0 +0.04250003832291116 0.8719799289625894 0 +0.08453617853796393 0.849490431102492 0 +0.4568591744667994 0.8749646885891479 0 +0.4091619142239268 0.9075571612900571 0 +0.04311519436129138 0.1253222927888206 0 +0.08638545191260635 0.1003760082534542 0 +0.08143951686858684 0.04605799348813314 0 +0.1299038105679871 0.07500000000040882 0 +0.1291547065148689 0.1238278498046576 0 +0.1729617141100337 0.09705533780903358 0 +0.2210475464026823 0.07311799749200407 0 +0.2121410437203431 0.1240404098500114 0 +0.2609214707917657 0.09783090151366602 0 +0.9577548213270598 0.3749900213278653 0 +0.1969149298890575 0.7249990465655037 0 +0.5249999999987085 0.04330127018919272 0 +0.5499999999985664 0.08660254037849818 0 +0.5999999999986141 0.08660254037874301 0 +0.5732076486424678 0.1285232051163829 0 +0.6271646611772891 0.133731112738492 0 +0.5197570975215425 0.1268419949038727 0 +0.6503607768618354 0.08724042407403498 0 +0.6805075989938938 0.1298619253837945 0 +0.7009780626418696 0.08670187346421093 0 +0.723706544204836 0.1279260804164832 0 +0.7560547114456557 0.09343760599329914 0 +0.8746806644031563 0.04307509652737851 0 +0.9566987298107054 0.2749999999994778 0 +0.6250601294760378 0.04340758413860749 0 +0.2401958109557048 0.3500059445558459 0 +0.2775135816834995 0.8773048153269645 0 +0.4566987298107762 0.5749999999999621 0 +0.413397459621575 0.599999999999927 0 +0.3700961894323919 0.5749999999998715 0 +0.3700961894324183 0.5249999999999589 0 +0.04127871425268339 0.3749249988936385 0 +0.4566987298107694 0.6249999999999646 0 +0.4133974596214428 0.6499999999997883 0 +0.2820313857775241 0.8247700375995279 0 +0.1969246305519431 0.6250001252370179 0 +0.08569757301204467 0.148686206903467 0 +0.1235154449242572 0.1733669783603554 0 +0.4566987298107781 0.7749999999999996 0 +0.3268592261749324 0.2001113828734097 0 +0.3701069072543669 0.1750185638126835 0 +0.4133992459252066 0.2000030939691757 0 +0.4133995436424987 0.1500036096307278 0 +0.3700983230729292 0.1250036955743397 0 +0.4127931361227511 0.1010491540728449 0 +0.4643847186991359 0.1256409747858049 0 +0.3709057972876633 0.07609655596763798 0 +0.3295361766483374 0.1012516567994244 0 +0.4177837134186551 0.04849267841999051 0 +0.7752626818440511 0.04332805082611701 0 +0.7262159093215194 0.04446146674553049 0 +0.3267949192432332 0.5499999999998454 0 +0.1970012858956514 0.5750001566211709 0 +0.456698729810777 0.4243554144234847 0 +0.4996298938899005 0.3980815905459545 0 +0.8031238207224871 0.2249745888704421 0 +0.6253222927882887 0.4568848056388577 0 +0.6003760082530059 0.4136145480876561 0 +0.6497757920799156 0.4136612942226305 0 +0.6250253000554828 0.3701763432768798 0 +0.6767122530443329 0.3714444895119191 0 +0.6528267423246743 0.3280260263682412 0 +0.7271668956210002 0.3766449766675489 0 +0.6015912417386127 0.3275694086808839 0 +0.546005121186425 0.4177339227153298 0 +0.575333476614867 0.3709520170446267 0 +0.6746189908132442 0.4571187906077035 0 +0.2834936490541751 0.4250000000001182 0 +0.7750643281474813 0.4595979569447049 0 +0.3279969479038325 0.1514992513288905 0 +0.270294436107 0.04584575151830302 0 +0.4133977573388259 0.2500005156618333 0 +0.8897114317026038 0.3250000000002788 0 +0.3267949192433037 0.4000000000001361 0 +0.04066270370297402 0.7764941934974348 0 +0.3700961894322892 0.6249999999997236 0 +0.3267949192432948 0.3500000000002306 0 +0.9566987298109041 0.1749999999995603 0 +0.4566987298107167 0.3750000000001421 0 +0.4565669311903869 0.2252551572761308 0 +0.3264342162508831 0.9494810250322635 0 +0.8039085986558747 0.1745700149423582 0 +0.1128170624400073 0.3236670214159245 0 +0.2402182048811268 0.3000447319726176 0 +0.3267949192430656 0.6999999999993651 0 +0.4978683978955004 0.08449594494604343 0 +0.4566987298106754 0.1750000000005258 0 +0.4989918945447412 0.3499381728500433 0 +0.4566987298107406 0.724999999999867 0 +0.3267949192431943 0.5999999999997128 0 +0.1968912828288999 0.4750003016420165 0 +0.4125787192934603 0.8510119220548549 0 +0.4132610029001055 0.8001686536753924 0 +0.2834942210692822 0.3750009907594891 0 +0.1536278364729008 0.5000003739640793 0 +0.240194258915888 0.4000032563438776 0 +0.1111926989368446 0.5250000268031971 0 +0.3700961894321958 0.6749999999995572 0 +0.1968930902032182 0.4250034321063578 0 +0.4133974596214268 0.699999999999642 0 +0.5448748726224287 0.1721664850222873 0 +0.7118235378800305 0.1716304707452299 0 +0.6997653361408653 0.4153578216014462 0 +0.6296565585981726 0.1818819788103112 0 +0.4254372001960948 0.9583812330199448 0 +0.07661977849406189 0.9596291756852154 0 +0.95307282801139 0.07795330506958331 0 +0.923373825580313 0.4554561683103004 0 +0.7488394277439033 0.4187318053084058 0 +0.7620497328870218 0.3461165485294553 0 +0.8999999999997218 0.08097344149834382 0 +0.7633265291612522 0.1493933528951607 0 +0.8013852629994044 0.0790431828483063 0 +0.8022377817800301 0.4185216713896819 0 +0.1936393345435158 0.8693677619870136 0 +0.4530685711562871 0.08128761922979176 0 +0.6765690845259915 0.2920710631917504 0 +0.2875015636282887 0.1319419116753489 0 +0.2483236425138556 0.1511807191013859 0 +0.8749999999999996 0.4652238251818912 0 +0.0390237343369057 0.6824089094277546 0 +0.07990472122704821 0.8005132877582968 0 +0.074843052028586 0.2988345477665329 0 +0.1667080418465754 0.1503489020178126 0 +0.1606893797559761 0.2005261540516154 0 +0.03902702860755038 0.6198272581762772 0 +0.08190980268676132 0.1939078376855576 0 +0.03905847369468828 0.2233439794436273 0 +0.3740857608536885 0.033917508872546 0 +0.07712153967911747 0.7501967849206319 0 +0.5334936490538766 0.3749999999999999 0 +0.1255665986823087 0.03339040557519336 0 +0.9232050807570187 0.1499999999996327 0 +0.9232050807570018 0.1999999999995609 0 +0.07679491924292095 0.5000000000018661 0 +0.07694589153052413 0.5499464508380731 0 +0.9232050807569379 0.2499999999995409 0 +0.9235419220874609 0.4003550923214798 0 +0.9232896855212954 0.3500209617135419 0 +0.9232191815508012 0.3000034936187781 0 +0.07679491924331235 0.4500000000012652 0 +0.07585892848436257 0.3494240925085729 0 +0.0766665635803451 0.399971760525902 0 +0.07687010113550455 0.5999999999998198 0 +0.1196039060330921 0.8250756563499276 0 +0.1148049474275764 0.8719857937560548 0 +0.07687010113628337 0.6999999999993193 0 +0.03812568041634121 0.2737420921931625 0 +0.3621762666320555 0.9168351706308571 0 +0.4949759526415395 0.1587019052846017 0 +0.224335414180262 0.03264637092141962 0 +0.2905208718860268 0.9188702367909223 0 +0.5811297632093186 0.2905208718852724 0 +0.2947679190005745 0.08011193437013248 0 +0.4633974596215362 0.9633974596215914 0 +0.03660254037856588 0.03660254037864849 0 +0.9633974596214215 0.03660254037851846 0 +0.03660254037858465 0.9633974596214154 0 +0.5366025403784426 0.4633974596215554 0 +0.9636520367722665 0.4636520367722302 0 +0.5955856378236344 0.1625047505404797 0 +0.7663805979810372 0.3835854935149701 0 +0.9251253471385955 0.1123767883649336 0 +0.2309867184120597 0.8850989416081687 0 +0.847230251594835 0.4328765465229614 0 +0.6642208584515532 0.1640425812304157 0 +0.07791370977919647 0.8849505039432525 0 +$EndNodes +$Elements +7 814 1 814 +1 1 1 20 +1 1 7 +2 7 8 +3 8 9 +4 9 10 +5 10 11 +6 11 12 +7 12 13 +8 13 14 +9 14 15 +10 15 16 +11 16 17 +12 17 18 +13 18 19 +14 19 20 +15 20 21 +16 21 22 +17 22 23 +18 23 24 +19 24 25 +20 25 2 +1 2 1 10 +21 2 26 +22 26 27 +23 27 28 +24 28 29 +25 29 30 +26 30 31 +27 31 32 +28 32 33 +29 33 34 +30 34 3 +1 3 1 10 +31 3 35 +32 35 36 +33 36 37 +34 37 38 +35 38 39 +36 39 40 +37 40 41 +38 41 42 +39 42 43 +40 43 4 +1 4 1 10 +41 4 44 +42 44 45 +43 45 46 +44 46 47 +45 47 48 +46 48 49 +47 49 50 +48 50 51 +49 51 52 +50 52 5 +1 5 1 10 +51 5 53 +52 53 54 +53 54 55 +54 55 56 +55 56 57 +56 57 58 +57 58 59 +58 59 60 +59 60 61 +60 61 6 +1 6 1 20 +61 6 62 +62 62 63 +63 63 64 +64 64 65 +65 65 66 +66 66 67 +67 67 68 +68 68 69 +69 69 70 +70 70 71 +71 71 72 +72 72 73 +73 73 74 +74 74 75 +75 75 76 +76 76 77 +77 77 78 +78 78 79 +79 79 80 +80 80 1 +2 1 2 734 +81 258 287 327 +82 168 302 303 +83 263 221 353 +84 149 241 346 +85 241 149 390 +86 287 258 391 +87 220 109 264 +88 15 85 290 +89 250 248 312 +90 225 178 356 +91 168 303 358 +92 362 68 367 +93 365 172 366 +94 279 365 366 +95 221 263 354 +96 235 362 367 +97 220 264 352 +98 262 263 353 +99 302 168 304 +100 312 248 392 +101 7 8 244 +102 14 15 290 +103 194 296 306 +104 145 303 305 +105 178 225 387 +106 260 262 343 +107 224 225 356 +108 290 85 357 +109 81 288 289 +110 245 112 247 +111 170 305 307 +112 250 359 360 +113 7 244 397 +114 194 306 400 +115 249 250 360 +116 260 343 407 +117 247 112 248 +118 109 220 354 +119 241 195 346 +120 322 206 390 +121 155 343 353 +122 145 305 394 +123 349 219 379 +124 183 228 347 +125 204 216 368 +126 368 216 369 +127 185 349 379 +128 258 256 342 +129 196 348 352 +130 212 220 352 +131 364 324 383 +132 363 316 371 +133 76 369 389 +134 228 183 408 +135 249 172 365 +136 327 287 357 +137 343 262 353 +138 247 248 249 +139 206 322 393 +140 305 170 394 +141 149 322 390 +142 303 145 358 +143 172 249 360 +144 352 348 404 +145 342 256 402 +146 302 304 344 +147 268 226 405 +148 249 248 250 +149 288 81 370 +150 170 307 372 +151 112 245 373 +152 182 268 405 +153 369 216 389 +154 82 364 383 +155 198 363 371 +156 304 168 351 +157 240 241 333 +158 241 205 333 +159 291 263 292 +160 263 261 292 +161 261 91 292 +162 108 240 333 +163 259 91 261 +164 90 330 341 +165 280 200 330 +166 280 108 334 +167 57 58 223 +168 94 57 223 +169 91 259 266 +170 108 333 334 +171 223 224 226 +172 200 280 334 +173 63 64 238 +174 64 136 238 +175 71 72 222 +176 94 223 226 +177 214 202 340 +178 236 157 294 +179 30 31 265 +180 253 142 254 +181 59 60 227 +182 18 19 266 +183 277 233 294 +184 95 30 265 +185 158 231 277 +186 70 71 237 +187 71 222 237 +188 58 124 223 +189 84 70 237 +190 62 63 183 +191 233 236 294 +192 330 200 341 +193 231 233 277 +194 267 190 325 +195 271 317 331 +196 229 158 252 +197 46 90 274 +198 270 275 317 +199 158 229 231 +200 175 267 325 +201 274 90 275 +202 242 184 243 +203 58 59 124 +204 72 134 222 +205 238 136 239 +206 142 18 266 +207 223 124 224 +208 270 274 275 +209 243 184 244 +210 79 184 242 +211 19 91 266 +212 255 142 266 +213 183 63 238 +214 259 255 266 +215 242 243 278 +216 17 142 253 +217 108 50 240 +218 125 32 251 +219 317 174 331 +220 48 280 330 +221 101 96 272 +222 243 246 278 +223 24 196 264 +224 243 244 245 +225 157 123 294 +226 213 214 215 +227 240 195 241 +228 159 187 252 +229 45 46 274 +230 182 147 276 +231 332 157 336 +232 268 182 276 +233 187 229 252 +234 96 92 272 +235 134 73 273 +236 31 125 265 +237 278 246 279 +238 157 236 336 +239 147 141 276 +240 51 195 240 +241 49 108 280 +242 79 80 184 +243 51 52 195 +244 272 271 293 +245 213 202 214 +246 201 190 267 +247 106 101 293 +248 159 177 187 +249 124 59 227 +250 187 177 198 +251 103 78 242 +252 189 205 206 +253 141 166 276 +254 31 32 125 +255 254 256 258 +256 22 109 291 +257 187 198 207 +258 246 245 247 +259 101 272 293 +260 33 185 251 +261 105 101 106 +262 24 25 196 +263 254 142 255 +264 44 45 269 +265 122 106 123 +266 109 23 264 +267 72 73 134 +268 254 255 256 +269 126 106 293 +270 191 190 201 +271 158 131 252 +272 133 159 252 +273 281 169 311 +274 140 113 281 +275 103 242 278 +276 21 291 292 +277 93 92 96 +278 122 123 157 +279 89 88 92 +280 232 230 235 +281 49 50 108 +282 159 171 177 +283 269 45 274 +284 113 138 281 +285 78 79 242 +286 105 106 122 +287 199 180 208 +288 164 270 271 +289 123 127 294 +290 140 281 282 +291 77 78 103 +292 33 34 185 +293 138 169 281 +294 93 98 295 +295 189 188 205 +296 293 271 331 +297 32 33 251 +298 93 96 97 +299 74 82 273 +300 89 92 93 +301 92 164 272 +302 92 88 164 +303 23 24 264 +304 129 131 158 +305 22 23 109 +306 164 271 272 +307 129 158 277 +308 48 49 280 +309 138 148 169 +310 50 51 240 +311 127 129 277 +312 89 93 295 +313 4 88 89 +314 229 207 230 +315 305 301 307 +316 131 133 252 +317 17 18 142 +318 166 188 189 +319 21 22 291 +320 73 74 273 +321 189 206 268 +322 100 101 105 +323 191 201 202 +324 180 160 297 +325 4 44 88 +326 46 47 90 +327 166 189 276 +328 127 277 294 +329 100 96 101 +330 190 162 325 +331 91 20 292 +332 332 214 340 +333 299 87 306 +334 209 180 297 +335 176 162 190 +336 179 160 180 +337 286 288 290 +338 123 126 127 +339 98 97 99 +340 208 180 209 +341 253 254 327 +342 286 285 288 +343 284 285 286 +344 187 207 229 +345 299 306 307 +346 298 87 299 +347 243 245 246 +348 284 282 285 +349 20 21 292 +350 231 230 232 +351 4 89 194 +352 102 99 104 +353 74 75 82 +354 123 106 126 +355 303 301 305 +356 171 147 182 +357 301 299 307 +358 164 269 270 +359 194 89 295 +360 42 87 298 +361 88 44 269 +362 107 104 110 +363 140 282 283 +364 42 43 87 +365 256 255 257 +366 127 128 129 +367 19 20 91 +368 179 180 199 +369 110 113 140 +370 120 118 163 +371 16 17 253 +372 208 209 210 +373 283 282 284 +374 224 124 225 +375 133 132 141 +376 177 171 178 +377 146 152 153 +378 153 152 154 +379 102 104 107 +380 300 299 301 +381 260 261 262 +382 296 295 320 +383 145 144 146 +384 302 301 303 +385 93 97 98 +386 41 298 308 +387 121 144 145 +388 98 99 102 +389 298 300 308 +390 160 154 297 +391 41 42 298 +392 141 132 165 +393 298 299 300 +394 40 41 308 +395 141 165 166 +396 176 190 191 +397 129 130 131 +398 163 118 321 +399 194 295 296 +400 110 140 313 +401 133 147 159 +402 300 301 302 +403 99 97 315 +404 153 154 156 +405 97 96 100 +406 288 285 289 +407 15 16 85 +408 86 40 308 +409 146 151 152 +410 105 122 137 +411 164 88 269 +412 173 181 193 +413 229 230 231 +414 282 281 311 +415 285 282 311 +416 233 232 234 +417 117 116 118 +418 154 155 297 +419 189 268 276 +420 262 261 263 +421 156 154 160 +422 107 110 313 +423 39 40 86 +424 317 275 339 +425 313 283 321 +426 127 126 128 +427 116 313 321 +428 37 38 135 +429 159 147 171 +430 154 152 155 +431 116 107 313 +432 12 13 81 +433 115 107 116 +434 38 39 310 +435 104 99 318 +436 129 128 130 +437 270 269 274 +438 140 283 313 +439 225 227 228 +440 166 186 188 +441 231 232 233 +442 133 141 147 +443 121 143 144 +444 284 286 287 +445 153 156 168 +446 193 181 217 +447 131 130 132 +448 257 259 260 +449 102 107 115 +450 135 38 310 +451 115 116 117 +452 271 270 317 +453 138 139 148 +454 203 202 213 +455 289 285 311 +456 11 12 312 +457 39 86 310 +458 146 144 151 +459 148 161 172 +460 156 160 173 +461 113 114 138 +462 99 315 318 +463 117 118 119 +464 114 111 175 +465 110 104 111 +466 131 132 133 +467 66 67 83 +468 173 160 179 +469 69 70 84 +470 297 155 323 +471 173 179 181 +472 110 111 113 +473 295 98 320 +474 217 218 219 +475 117 119 170 +476 29 30 95 +477 105 137 309 +478 97 100 315 +479 65 66 316 +480 56 57 94 +481 12 81 312 +482 233 234 236 +483 64 65 136 +484 283 284 328 +485 217 181 218 +486 321 283 328 +487 179 199 314 +488 214 332 336 +489 210 211 212 +490 119 118 120 +491 119 120 121 +492 260 259 261 +493 98 102 320 +494 209 297 323 +495 175 111 318 +496 47 48 330 +497 100 105 309 +498 257 255 259 +499 28 29 319 +500 210 209 211 +501 130 128 174 +502 118 116 321 +503 166 165 186 +504 181 179 314 +505 218 181 314 +506 111 104 318 +507 121 120 143 +508 167 197 204 +509 296 320 329 +510 27 28 150 +511 163 321 328 +512 102 115 320 +513 55 56 322 +514 220 211 221 +515 186 165 192 +516 126 293 331 +517 132 130 326 +518 167 162 176 +519 191 202 203 +520 113 111 114 +521 138 114 139 +522 66 83 316 +523 165 132 326 +524 254 258 327 +525 100 309 315 +526 85 253 327 +527 275 90 341 +528 9 10 112 +529 215 214 336 +530 186 200 334 +531 148 139 161 +532 54 55 149 +533 114 175 325 +534 136 65 316 +535 186 192 200 +536 204 197 216 +537 117 170 329 +538 161 139 162 +539 212 211 220 +540 161 162 167 +541 130 174 326 +542 29 95 319 +543 167 176 197 +544 205 188 333 +545 192 165 326 +546 128 126 331 +547 188 186 334 +548 176 191 324 +549 137 122 332 +550 90 47 330 +551 56 94 322 +552 122 157 332 +553 115 117 329 +554 284 287 328 +555 320 115 329 +556 150 28 319 +557 211 209 323 +558 197 176 324 +559 162 139 325 +560 221 211 323 +561 225 124 227 +562 191 203 324 +563 267 175 335 +564 174 128 331 +565 139 114 325 +566 85 16 253 +567 149 55 322 +568 315 309 335 +569 335 309 337 +570 201 267 337 +571 202 201 340 +572 333 188 334 +573 318 315 335 +574 175 318 335 +575 267 335 337 +576 201 337 340 +577 174 317 339 +578 215 336 338 +579 309 137 337 +580 339 275 341 +581 326 174 339 +582 337 137 340 +583 236 234 338 +584 336 236 338 +585 137 332 340 +586 192 326 339 +587 192 339 341 +588 200 192 341 +589 355 217 406 +590 120 163 342 +591 86 344 350 +592 279 204 368 +593 204 279 366 +594 143 120 342 +595 289 359 395 +596 81 289 395 +597 155 152 343 +598 86 308 344 +599 300 302 344 +600 144 143 345 +601 60 61 347 +602 228 227 347 +603 53 54 346 +604 26 27 348 +605 247 249 365 +606 135 355 406 +607 344 304 350 +608 173 193 351 +609 168 156 351 +610 250 312 395 +611 152 151 343 +612 308 300 344 +613 323 155 353 +614 248 112 392 +615 263 291 354 +616 220 221 354 +617 193 217 355 +618 151 144 345 +619 171 182 356 +620 227 60 347 +621 54 149 346 +622 27 150 348 +623 286 290 357 +624 143 342 402 +625 145 146 358 +626 153 168 358 +627 289 311 359 +628 342 163 391 +629 148 172 360 +630 37 135 361 +631 35 36 349 +632 67 68 362 +633 136 316 363 +634 278 279 368 +635 197 324 364 +636 279 246 365 +637 68 69 367 +638 167 204 366 +639 76 77 369 +640 172 161 366 +641 81 13 370 +642 329 170 372 +643 9 112 373 +644 232 235 385 +645 198 177 386 +646 235 230 388 +647 82 75 389 +648 258 342 391 +649 205 241 390 +650 328 287 391 +651 11 312 392 +652 112 10 392 +653 359 250 395 +654 121 145 394 +655 170 119 394 +656 322 94 393 +657 219 349 361 +658 257 260 407 +659 5 53 396 +660 80 1 397 +661 1 7 397 +662 2 26 398 +663 25 2 398 +664 244 184 397 +665 52 5 396 +666 61 6 399 +667 6 62 399 +668 306 87 400 +669 4 194 400 +670 43 4 400 +671 310 86 350 +672 183 238 408 +673 156 173 351 +674 34 3 401 +675 3 35 401 +676 310 350 355 +677 85 327 357 +678 345 257 407 +679 349 36 361 +680 225 228 387 +681 264 196 352 +682 221 323 353 +683 291 109 354 +684 135 310 355 +685 178 171 356 +686 343 151 407 +687 355 350 403 +688 287 286 357 +689 14 290 370 +690 244 8 373 +691 306 296 372 +692 226 268 393 +693 146 153 358 +694 311 169 359 +695 362 235 388 +696 363 198 386 +697 235 367 385 +698 364 82 389 +699 239 386 387 +700 386 178 387 +701 359 169 360 +702 169 148 360 +703 36 37 361 +704 83 67 362 +705 239 136 363 +706 216 197 364 +707 246 247 365 +708 69 84 367 +709 161 167 366 +710 103 368 369 +711 103 278 368 +712 77 103 369 +713 312 81 395 +714 268 206 393 +715 371 83 388 +716 84 377 385 +717 213 382 384 +718 382 134 384 +719 380 314 381 +720 376 338 377 +721 125 380 381 +722 207 198 371 +723 338 234 377 +724 290 288 370 +725 199 378 381 +726 213 215 382 +727 199 208 378 +728 208 210 375 +729 210 212 374 +730 215 376 382 +731 8 9 373 +732 125 251 380 +733 215 338 376 +734 210 374 375 +735 13 14 370 +736 208 375 378 +737 134 273 384 +738 314 199 381 +739 251 185 379 +740 273 82 383 +741 251 379 380 +742 237 376 377 +743 316 83 371 +744 378 265 381 +745 273 383 384 +746 245 244 373 +747 95 265 378 +748 307 306 372 +749 84 237 377 +750 296 329 372 +751 237 222 376 +752 319 95 375 +753 150 319 374 +754 222 134 382 +755 374 319 375 +756 265 125 381 +757 219 218 379 +758 379 218 380 +759 218 314 380 +760 375 95 378 +761 376 222 382 +762 203 213 384 +763 324 203 383 +764 383 203 384 +765 207 371 388 +766 377 234 385 +767 234 232 385 +768 177 178 386 +769 230 207 388 +770 75 76 389 +771 206 205 390 +772 163 328 391 +773 256 257 402 +774 10 11 392 +775 94 226 393 +776 119 121 394 +777 193 355 403 +778 348 150 404 +779 347 61 399 +780 53 346 396 +781 26 348 398 +782 185 34 401 +783 195 52 396 +784 184 80 397 +785 196 25 398 +786 62 183 399 +787 87 43 400 +788 361 135 406 +789 345 143 402 +790 257 345 402 +791 226 224 405 +792 351 193 403 +793 356 182 405 +794 217 219 406 +795 238 239 408 +796 183 347 399 +797 346 195 396 +798 348 196 398 +799 349 185 401 +800 35 349 401 +801 83 362 388 +802 239 363 386 +803 216 364 389 +804 367 84 385 +805 350 304 403 +806 151 345 407 +807 304 351 403 +808 239 387 408 +809 387 228 408 +810 224 356 405 +811 212 352 404 +812 150 374 404 +813 374 212 404 +814 219 361 406 +$EndElements diff --git a/src/argyris_biharmonic.jl b/src/argyris_biharmonic.jl new file mode 100644 index 0000000..4a4d8c6 --- /dev/null +++ b/src/argyris_biharmonic.jl @@ -0,0 +1,125 @@ +# In this tutorial, we will learn +# - How to use the Argyris element, the classical $C^1$-conforming triangle +# - How to impose clamped boundary conditions when the degrees of freedom are derivatives +# - How to measure the $C^1$ continuity of a discrete space +# +# ## Problem statement +# +# We solve the biharmonic problem +# +# ```math +# \Delta^2 u = f \ \text{ in } \Omega = (0,1)^2, +# ``` +# +# with a manufactured solution, using an $H^2$-*conforming* discretisation. The +# variational form is +# +# ```math +# a(u,v) = \int_\Omega D^2 u : D^2 v, \qquad \ell(v) = \int_\Omega f\, v, +# ``` +# +# and, unlike in the Morley tutorial, here $D^2 u$ really is the distributional +# Hessian: the discrete space is a subspace of $H^2(\Omega)$. +# +# ## The element +# +# The **Argyris** triangle is the quintic space $P_5(K)$, of dimension 21, with +# per vertex the value, the two first derivatives and the three second +# derivatives, and per edge one normal-derivative moment: +# +# ```math +# u(v),\quad \partial_i u(v),\quad \partial_{ij} u(v), \qquad \int_e \nabla u\cdot n . +# ``` +# +# Sharing all six vertex quantities between the elements meeting at a vertex, and +# the normal-derivative moment between the two elements sharing an edge, makes the +# assembled space $C^1$. It is the cheapest $C^1$ triangle there is. +# +# Two consequences worth noting. The vertex degrees of freedom are *point +# evaluations of derivatives*, so they are not preserved by the pullback and the +# element needs a cell-dependent change of basis — Gridap does this internally. +# And because those degrees of freedom include second derivatives, a boundary +# condition on `dirichlet_tags` constrains the boundary Hessian too, which is why +# below we interpolate the exact solution rather than setting zero. + +using Gridap + +# ## Discrete model +# +# Argyris is defined on triangles, so we simplexify a Cartesian mesh. + +n = 8 +model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) + +# ## Manufactured solution +# +# We take $u = p(x)p(y)$ with $p(t) = t^2(1-t)^2$, which vanishes together with +# its gradient on $\partial\Omega$, and compute $f = \Delta^2 u$ by hand. + +p(t) = t^2 * (1 - t)^2 +pdd(t) = 2 - 12*t + 12*t^2 + +uex(x) = p(x[1]) * p(x[2]) +f(x) = 24*p(x[2]) + 2*pdd(x[1]) * pdd(x[2]) + 24*p(x[1]) + +# ## FE space +# +# `argyris` is a scalar element of order 5 on triangles. The boundary degrees of +# freedom include the vertex Hessians, and $D^2 u \neq 0$ on $\partial\Omega$ +# even for a clamped plate, so the Dirichlet data is the exact solution: its +# interpolation supplies the right value for every boundary degree of freedom at +# once. + +V = FESpace(model, ReferenceFE(argyris, Float64, 5); dirichlet_tags="boundary") +U = TrialFESpace(V, uex) + +Ω = Triangulation(model) +dΩ = Measure(Ω, 12) + +a(u, v) = ∫( ∇∇(u) ⊙ ∇∇(v) )dΩ +l(v) = ∫( f * v )dΩ + +op = AffineFEOperator(a, l, U, V) +uh = solve(op) + +# ## Errors +# +# The space is $P_5$, so we expect $O(h^6)$ in $L^2$ and $O(h^4)$ in the $H^2$ +# seminorm. + +Hex = ∇∇(uex) + +e = uh - uex +eH = ∇∇(uh) - Hex +l2 = sqrt(sum( ∫( e * e )dΩ )) +h2 = sqrt(sum( ∫( eH ⊙ eH )dΩ )) +println("L2 error: ", l2) +println("H2 seminorm error: ", h2) + +# ## The space really is $C^1$ +# +# Both the value and the *full* gradient are continuous across every interior +# edge — not just the mean normal derivative, as in Morley. We check on a +# function that is not in the space, so that the jumps are not zero for the +# trivial reason. Note the space here carries no `dirichlet_tags`: interpolating +# into a constrained space would zero the boundary degrees of freedom and inflate +# the jumps near $\partial\Omega$. + +w(x) = sin(2*x[1]) * cos(3*x[2]) +wh = interpolate(w, FESpace(model, ReferenceFE(argyris, Float64, 5))) + +Λ = SkeletonTriangulation(model) +dΛ = Measure(Λ, 12) + +println("value jump: ", sqrt(sum( ∫( jump(wh) * jump(wh) )dΛ ))) +println("gradient jump: ", sqrt(sum( ∫( jump(∇(wh)) ⋅ jump(∇(wh)) )dΛ ))) + +# ... but not $C^2$: the Hessian does jump, as it must for a finite element space +# built on a triangulation. + +println("Hessian jump: ", sqrt(sum( ∫( jump(∇∇(wh)) ⊙ jump(∇∇(wh)) )dΛ ))) + +# ## Visualisation + +mkpath("output_path") +writevtk(Ω, "output_path/argyris", cellfields=["u" => uh, "error" => e]) diff --git a/src/arnold_winther_elasticity.jl b/src/arnold_winther_elasticity.jl new file mode 100644 index 0000000..11ce70e --- /dev/null +++ b/src/arnold_winther_elasticity.jl @@ -0,0 +1,158 @@ +# In this tutorial, we will learn +# - How to solve linear elasticity in *mixed* (Hellinger--Reissner) form +# - How to use the Arnold--Winther elements, which impose the symmetry of the +# stress tensor exactly +# - What the conforming and nonconforming variants buy you +# +# ## Problem statement +# +# The Hellinger--Reissner formulation of linear elasticity takes the stress +# $\sigma$ and the displacement $u$ as unknowns: with the identity compliance, +# +# ```math +# \int_\Omega \sigma : \tau + b(\tau, u) = 0, \qquad +# b(\sigma, v) = -\int_\Omega f\cdot v, \qquad +# b(\tau, v) = \sum_K \int_K \operatorname{div}\tau \cdot v . +# ``` +# +# Its difficulty is well known: a stable pair needs a stress space that is +# **symmetric**, **$H(\operatorname{div})$-conforming** and matched to the +# displacement space. Elements satisfying all three eluded the field for decades. +# +# ## The elements +# +# Arnold and Winther gave the first two on triangles. +# +# The **conforming** element, `aw_c`, has 24 degrees of freedom: +# +# ```math +# AW_c(K) = \{\tau \in P_3(K;\mathbb{S}) : \operatorname{div}\tau \in P_1(K;\mathbb{R}^2)\}, +# ``` +# +# with the three components of $\tau$ at each vertex, the degree 0 and 1 moments +# of $n\cdot\tau n$ and $n\cdot\tau t$ on each edge, and three interior moments. +# It is genuinely $H(\operatorname{div};\mathbb{S})$-conforming: the full traction +# $\tau\cdot n$ is continuous. +# +# The **nonconforming** element, `aw_nc`, has 15: +# +# ```math +# AW_{nc}(K) = \{\tau \in P_2(K;\mathbb{S}) : (n\cdot\tau n)|_e \in P_1(e)\}, +# ``` +# +# with only the edge moments and three interior ones. The traction jumps, but the +# jump is orthogonal to the $P_1$ displacement space, which is enough for +# convergence — the standard trade of accuracy for cost. +# +# Both are *constrained* subspaces of a larger polynomial space, so for them +# `length(get_prebasis(reffe))` is larger than `num_dofs(reffe)`. That is by +# design and nothing downstream depends on the two being equal. + +using Gridap +using Gridap.TensorValues +using Gridap.MultiField + +# ## Discrete model + +n = 8 +model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) + +Ω = Triangulation(model) +Λ = SkeletonTriangulation(model) +Γ = BoundaryTriangulation(model) +degree = 8 +dΩ = Measure(Ω, degree) +dΛ = Measure(Λ, degree) +dΓ = Measure(Γ, degree) +nΛ = get_normal_vector(Λ) +nΓ = get_normal_vector(Γ) + +# ## Manufactured solution +# +# A displacement vanishing on $\partial\Omega$, its symmetric gradient, and the +# body force that produces it. + +const π2 = pi^2 +uex(x) = VectorValue(sin(pi*x[1]) * sin(pi*x[2]), sin(2*pi*x[1]) * sin(pi*x[2])) + +function σex(x) + u1x = pi * cos(pi*x[1]) * sin(pi*x[2]) + u1y = pi * sin(pi*x[1]) * cos(pi*x[2]) + u2x = 2*pi * cos(2*pi*x[1]) * sin(pi*x[2]) + u2y = pi * sin(2*pi*x[1]) * cos(pi*x[2]) + SymTensorValue{2,Float64}(u1x, (u1y + u2x) / 2, u2y) +end + +function fex(x) + u1 = sin(pi*x[1]) * sin(pi*x[2]) + u2 = sin(2*pi*x[1]) * sin(pi*x[2]) + d11u1, d22u1 = -π2*u1, -π2*u1 + d12u1 = π2 * cos(pi*x[1]) * cos(pi*x[2]) + d11u2, d22u2 = -4*π2*u2, -π2*u2 + d12u2 = 2*π2 * cos(2*pi*x[1]) * cos(pi*x[2]) + VectorValue(-(d11u1 + (d22u1 + d12u2)/2), -((d12u1 + d11u2)/2 + d22u2)) +end + +# ## The solver +# +# The displacement space is discontinuous vector $P_1$. The bilinear form $b$ is +# integrated by parts element by element, so no derivative of $\tau$ is ever +# taken and each side of a facet contributes with its own traction — which is +# what keeps the form consistent for the *nonconforming* stress space. The +# displacement boundary condition $u = 0$ is natural here, so no +# `dirichlet_tags` appear. + +function solve_elasticity(reffe_σ) + Vσ = FESpace(model, reffe_σ) + Vu = FESpace(model, ReferenceFE(lagrangian, VectorValue{2,Float64}, 1); conformity=:L2) + X = MultiFieldFESpace([Vσ, Vu]) + Y = MultiFieldFESpace([Vσ, Vu]) + + b(τ, v) = ∫( -(τ ⊙ ε(v)) )dΩ + + ∫( ((τ.⁺ ⋅ nΛ.⁺) ⋅ v.⁺) + ((τ.⁻ ⋅ nΛ.⁻) ⋅ v.⁻) )dΛ + + ∫( ((τ ⋅ nΓ) ⋅ v) )dΓ + A((σ, u), (τ, v)) = ∫( σ ⊙ τ )dΩ + b(τ, u) + b(σ, v) + L((τ, v)) = ∫( -(fex ⋅ v) )dΩ + + solve(AffineFEOperator(A, L, X, Y)) +end + +σc, uc = solve_elasticity(ReferenceFE(TRI, aw_c, Float64)) +σn, un = solve_elasticity(ReferenceFE(TRI, aw_nc, Float64)) + +# ## Errors +# +# The conforming element gives $O(h^3)$ in the stress and $O(h^2)$ in the +# displacement; the nonconforming one, being built on $P_2$ rather than $P_3$, +# gives $O(h)$ and $O(h^2)$. + +for (name, σh, uh) in (("conforming ", σc, uc), ("nonconforming", σn, un)) + eσ = σh - σex + eu = uh - uex + println(name, " |sigma| = ", sqrt(sum( ∫( eσ ⊙ eσ )dΩ )), + " |u| = ", sqrt(sum( ∫( eu ⋅ eu )dΩ ))) +end + +# ## Conforming or not +# +# The difference between the two spaces, made visible: interpolate a stress field +# that is in neither, and measure the jump of the traction $\tau\cdot n$ across +# the interior edges. For the conforming element it vanishes; for the +# nonconforming one it does not, though its normal--normal part still does. + +g(x) = SymTensorValue{2,Float64}(sin(2*x[1]), cos(3*x[2]), sin(x[1] + x[2])) +nplus = get_normal_vector(Λ).⁺ + +for (name, reffe) in (("conforming ", ReferenceFE(TRI, aw_c, Float64)), + ("nonconforming", ReferenceFE(TRI, aw_nc, Float64))) + gh = interpolate(g, FESpace(model, reffe)) + jt = sqrt(sum( ∫( (jump(gh) ⋅ nplus) ⋅ (jump(gh) ⋅ nplus) )dΛ )) + jnn = sqrt(sum( ∫( (nplus ⋅ jump(gh) ⋅ nplus) * (nplus ⋅ jump(gh) ⋅ nplus) )dΛ )) + println(name, " traction jump = ", jt, " nn jump = ", jnn) +end + +# ## Visualisation + +mkpath("output_path") +writevtk(Ω, "output_path/arnold_winther", + cellfields=["sigma_c" => σc, "u_c" => uc, "sigma_nc" => σn, "u_nc" => un]) diff --git a/src/gls_stokes.jl b/src/gls_stokes.jl new file mode 100644 index 0000000..52f4bb8 --- /dev/null +++ b/src/gls_stokes.jl @@ -0,0 +1,130 @@ +# In this tutorial, we will learn +# - How to use the Gopalakrishnan--Lederer--Schöberl element, a +# *traceless*-matrix-valued element that is normal--tangential continuous +# - How to reconstruct a velocity gradient in it +# - Where it fits in the mass-conserving mixed stress formulation of Stokes +# +# ## The element +# +# The **GLS** element (of the second kind) is the space of *traceless* matrix +# fields +# +# ```math +# GLS_r(K) = P_r(K;\mathbb{M}_0), \qquad \mathbb{M}_0 = \{M : \operatorname{tr} M = 0\}, +# ``` +# +# of dimension $3(r+1)(r+2)/2$ for any $r \ge 0$, with degrees of freedom +# +# ```math +# \int_e (t\cdot M n)\,\mu_i \ \text{ on each edge}, \qquad +# \int_K (t_e\cdot M n_e)\, q \ \text{ in the cell}, +# ``` +# +# both of them the *same* normal--tangential functional — over the edges, which +# glue, and over the cell, which does not. Only $t\cdot M n$ is shared between two +# triangles. +# +# It completes a small family. Three elements share exactly the same shape of +# change of basis, because each pairs its matrix with two directions carried +# dually by its own push-forward: +# +# | element | degree of freedom | continuity | +# |---|---|---| +# | `hhj` | $n\cdot M n$ | normal--normal | +# | `regge` | $t\cdot M t$ | tangential--tangential | +# | `gls` | $t\cdot M n$ | normal--tangential | +# +# ## Why traceless +# +# In the mass-conserving mixed stress (MCS) formulation of Stokes flow, the +# unknowns are the velocity $u$ in an $H(\operatorname{div})$-conforming space, +# the pressure, and the *deviatoric* velocity gradient +# $\sigma = \nabla u - \tfrac{1}{2}(\operatorname{div} u) I$. Since +# $\operatorname{div} u = 0$ is imposed exactly by the velocity space, $\sigma$ is +# traceless — and only its normal--tangential component needs to be +# single-valued, because that is all the formulation tests. GLS is the space that +# supplies exactly that and no more. +# +# Note the element is built as a *constrained* subspace of the full matrix space +# $P_r(K;\mathbb{M})$, so `length(get_prebasis(reffe))` exceeds `num_dofs(reffe)` +# by the number of trace constraints. + +using Gridap +using Gridap.TensorValues +using Gridap.ReferenceFEs + +# ## Discrete model and space + +n = 8 +model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) + +Ω = Triangulation(model) +dΩ = Measure(Ω, 8) + +r = 1 +V = FESpace(model, ReferenceFE(gls, Float64, r)) + +reffe = ReferenceFE(TRI, gls, Float64, r) +println("dofs per cell: ", num_dofs(reffe), " ambient prebasis: ", length(get_prebasis(reffe))) + +# ## Reconstructing a deviatoric gradient +# +# Take a divergence-free velocity and project its gradient into the GLS space. +# Because $\operatorname{div} u = 0$, that gradient is already traceless, so it +# lies in the space as soon as it is a polynomial of degree $\le r$ — and then +# the projection is exact. Here $u = (y^2, -x^2)$, whose gradient is linear. + +u(x) = VectorValue(x[2]^2, -x[1]^2) +G = ∇(u) + +a(σ, τ) = ∫( σ ⊙ τ )dΩ +l(τ) = ∫( G ⊙ τ )dΩ +σh = solve(AffineFEOperator(a, l, V, V)) + +e = σh - G +println("L2 error of the projected gradient: ", sqrt(sum( ∫( e ⊙ e )dΩ ))) + +# For a velocity whose gradient is *not* in the space, the projection is the best +# approximation in $L^2$ rather than the field itself: + +ψ(x) = (x[1]^2 * (1 - x[1])^2) * (x[2]^2 * (1 - x[2])^2) +v(x) = VectorValue(∇(ψ)(x)[2], -∇(ψ)(x)[1]) +Gv = ∇(v) +σv = solve(AffineFEOperator(a, τ -> ∫( Gv ⊙ τ )dΩ, V, V)) +ev = σv - Gv +println("L2 error, gradient outside the space: ", sqrt(sum( ∫( ev ⊙ ev )dΩ ))) + +# The reconstruction really is traceless, pointwise: + +tr_σ = Operation(tr)(σh) +println("L2 norm of the trace: ", sqrt(sum( ∫( tr_σ * tr_σ )dΩ ))) + +# ## Normal--tangential continuity, and nothing more +# +# The defining property. As before we interpolate a field that is not in the +# space, so the full jump is not zero for a trivial reason. + +g(x) = TensorValue(sin(2*x[1]), cos(3*x[2]), sin(x[1] + x[2]), -sin(2*x[1])) +gh = interpolate(g, V) + +Λ = SkeletonTriangulation(model) +dΛ = Measure(Λ, 8) +nΛ = get_normal_vector(Λ).⁺ +tΛ = Operation(n -> VectorValue(n[2], -n[1]))(nΛ) + +nt = sqrt(sum( ∫( (tΛ ⋅ jump(gh) ⋅ nΛ) * (tΛ ⋅ jump(gh) ⋅ nΛ) )dΛ )) +full = sqrt(sum( ∫( jump(gh) ⊙ jump(gh) )dΛ )) +println("jump of the nt component: ", nt) +println("jump of the full matrix: ", full) + +# ## The family + +for rr in 0:2 + Vr = FESpace(model, ReferenceFE(gls, Float64, rr)) + println("r = ", rr, " dofs = ", num_free_dofs(Vr)) +end + +# ## Visualisation + +mkpath("output_path") +writevtk(Ω, "output_path/gls", cellfields=["sigma" => σh]) diff --git a/src/hhj_plate.jl b/src/hhj_plate.jl new file mode 100644 index 0000000..bde65c8 --- /dev/null +++ b/src/hhj_plate.jl @@ -0,0 +1,166 @@ +# In this tutorial, we will learn +# - How to solve the Kirchhoff plate with a *mixed* method +# - How to use the Hellan--Herrmann--Johnson element, a symmetric-tensor-valued +# element that is only normal--normal continuous +# - How to write a form that couples a cell integral to a skeleton integral +# +# ## Problem statement +# +# The clamped Kirchhoff plate again — $\Delta^2 w = f$ with $w = \partial_n w = 0$ +# — but written as a first-order system in the bending moment tensor +# $\sigma = D^2 w$: +# +# ```math +# \sigma - D^2 w = 0, \qquad -\operatorname{div}\operatorname{div}\sigma = f . +# ``` +# +# The advantage over the primal form is that the displacement space no longer +# needs to be $C^1$: continuous $P_{r+1}$ Lagrange suffices. +# +# ## The element +# +# The **Hellan--Herrmann--Johnson** element is the full symmetric-matrix-valued +# polynomial space $P_r(K;\mathbb{S})$, of any degree $r \ge 0$, with degrees of +# freedom +# +# ```math +# \int_e (n\cdot\sigma n)\,\mu_i \ \text{ on each edge}, \qquad +# \int_K \sigma : \tau \ \text{ in the cell}. +# ``` +# +# Only the **normal--normal** component $n\cdot\sigma n$ is shared between two +# elements; the rest of the tensor jumps. That is exactly the regularity the +# mixed plate formulation needs, and no more — the trade-off that makes the +# method cheap. +# +# Because $n\cdot\sigma n$ is quadratic in $n$, the element needs no normal sign +# convention at all: reversing an edge leaves the functional unchanged. +# +# ## The mixed form +# +# Following Arnold and Walker, with $W_h$ the continuous Lagrange space of degree +# $r+1$, +# +# ```math +# b(\varphi, v) = -\sum_K \int_K \varphi : D^2 v +# + \sum_e \int_e \varphi_{nn}\,[\![\partial_n v]\!], +# ``` +# +# the discrete problem is: find $(\sigma_h, w_h)$ with +# $(\sigma_h,\tau) + b(\tau,w_h) = 0$ and $b(\sigma_h,v) = -\langle f,v\rangle$. +# Only $\varphi_{nn}$ appears in the skeleton term, which is precisely the +# component HHJ makes single-valued. + +using Gridap +using Gridap.TensorValues +using Gridap.MultiField +using GridapGmsh + +# ## Discrete model and spaces + +n = 8 +model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) +r = 1 + +Ω = Triangulation(model) +Λ = SkeletonTriangulation(model) +Γ = BoundaryTriangulation(model) + +degree = 2*r + 4 +dΩ = Measure(Ω, degree) +dΛ = Measure(Λ, degree) +dΓ = Measure(Γ, degree) +nΛ = get_normal_vector(Λ).⁺ +nΓ = get_normal_vector(Γ) + +# `hhj` is a symmetric-tensor-valued element of degree `r`; the displacement uses +# ordinary Lagrange of degree `r+1`, with $w = 0$ imposed strongly. The other +# clamped condition, $\partial_n w = 0$, is imposed *weakly*, by including the +# boundary edges in the skeleton term with $[\![\partial_n v]\!] = \partial_n v$. + +Vσ = FESpace(model, ReferenceFE(hhj, Float64, r)) +Vw = FESpace(model, ReferenceFE(lagrangian, Float64, r + 1); dirichlet_tags="boundary") + +X = MultiFieldFESpace([Vσ, TrialFESpace(Vw, 0.0)]) +Y = MultiFieldFESpace([Vσ, Vw]) + +# ## Manufactured solution + +p(t) = t^2 * (1 - t)^2 +pdd(t) = 2 - 12*t + 12*t^2 + +wex(x) = p(x[1]) * p(x[2]) +σex = ∇∇(wex) +f(x) = 24*p(x[2]) + 2*pdd(x[1]) * pdd(x[2]) + 24*p(x[1]) + +# ## Weak form and solve + +bh(φ, v) = ∫( -(φ ⊙ ∇∇(v)) )dΩ + + ∫( (nΛ ⋅ mean(φ) ⋅ nΛ) * (jump(∇(v)) ⋅ nΛ) )dΛ + + ∫( (nΓ ⋅ φ ⋅ nΓ) * (∇(v) ⋅ nΓ) )dΓ + +A((σ, w), (τ, v)) = ∫( σ ⊙ τ )dΩ + bh(τ, w) + bh(σ, v) +L((τ, v)) = ∫( -(f * v) )dΩ + +σh, wh = solve(AffineFEOperator(A, L, X, Y)) + +# ## Errors +# +# The method converges at $O(h^{r+2})$ in the displacement and $O(h^{r+1})$ in +# the moment tensor. + +ew = wh - wex +eσ = σh - σex +println("L2 error, displacement: ", sqrt(sum( ∫( ew * ew )dΩ ))) +println("L2 error, moment: ", sqrt(sum( ∫( eσ ⊙ eσ )dΩ ))) + +# ## Normal--normal continuity, and nothing more +# +# The defining property of the space. We interpolate a tensor field that is not +# in it, so that the full jump is not zero for a trivial reason: the +# normal--normal component is continuous to round-off, while the tensor itself +# jumps at $O(1)$. + +g(x) = SymTensorValue{2,Float64}(sin(2*x[1]), cos(3*x[2]), sin(x[1] + x[2])) +gh = interpolate(g, Vσ) + +nn = sqrt(sum( ∫( (nΛ ⋅ jump(gh) ⋅ nΛ) * (nΛ ⋅ jump(gh) ⋅ nΛ) )dΛ )) +full = sqrt(sum( ∫( jump(gh) ⊙ jump(gh) )dΛ )) +println("jump of the nn component: ", nn) +println("jump of the full tensor: ", full) + +# ## A plate with a re-entrant corner +# +# Nothing above used the structure of the square. The same driver runs on an +# unstructured mesh of the L-shaped domain, where the clamped plate has a genuine +# corner singularity and no closed-form solution — a uniform load is enough to +# see it. + +model_L = GmshDiscreteModel("../models/lshape.msh") + +Ω_L = Triangulation(model_L) +Λ_L = SkeletonTriangulation(model_L) +Γ_L = BoundaryTriangulation(model_L) +dΩ_L = Measure(Ω_L, degree) +dΛ_L = Measure(Λ_L, degree) +dΓ_L = Measure(Γ_L, degree) +nΛ_L = get_normal_vector(Λ_L).⁺ +nΓ_L = get_normal_vector(Γ_L) + +Vσ_L = FESpace(model_L, ReferenceFE(hhj, Float64, r)) +Vw_L = FESpace(model_L, ReferenceFE(lagrangian, Float64, r + 1); dirichlet_tags="boundary") +X_L = MultiFieldFESpace([Vσ_L, TrialFESpace(Vw_L, 0.0)]) +Y_L = MultiFieldFESpace([Vσ_L, Vw_L]) + +bh_L(φ, v) = ∫( -(φ ⊙ ∇∇(v)) )dΩ_L + + ∫( (nΛ_L ⋅ mean(φ) ⋅ nΛ_L) * (jump(∇(v)) ⋅ nΛ_L) )dΛ_L + + ∫( (nΓ_L ⋅ φ ⋅ nΓ_L) * (∇(v) ⋅ nΓ_L) )dΓ_L + +A_L((σ, w), (τ, v)) = ∫( σ ⊙ τ )dΩ_L + bh_L(τ, w) + bh_L(σ, v) +L_L((τ, v)) = ∫( -(1.0 * v) )dΩ_L + +σh_L, wh_L = solve(AffineFEOperator(A_L, L_L, X_L, Y_L)) + +mkpath("output_path") +writevtk(Ω, "output_path/hhj_square", cellfields=["w" => wh, "sigma" => σh]) +writevtk(Ω_L, "output_path/hhj_lshape", cellfields=["w" => wh_L, "sigma" => σh_L]) diff --git a/src/morley_biharmonic.jl b/src/morley_biharmonic.jl new file mode 100644 index 0000000..e3db8fc --- /dev/null +++ b/src/morley_biharmonic.jl @@ -0,0 +1,128 @@ +# In this tutorial, we will learn +# - How to solve a fourth-order (biharmonic) problem in Gridap +# - How to use the Morley element, a *nonconforming* plate element +# - How to read an unstructured triangular mesh generated with Gmsh +# +# ## Problem statement +# +# We solve the **clamped plate**: find the transverse deflection $w$ of a thin +# plate $\Omega$ loaded by $f$ and clamped along its whole boundary, +# +# ```math +# \Delta^2 w = f \ \text{ in } \Omega, \qquad +# w = \frac{\partial w}{\partial n} = 0 \ \text{ on } \partial\Omega . +# ``` +# +# We take $\Omega$ the unit disk and $f \equiv 1$, for which the deflection is +# radially symmetric and known in closed form, +# +# ```math +# w(r) = \frac{(1-r^2)^2}{64}, +# ``` +# so the centre deflection is exactly $1/64$. +# +# ## Why a special element +# +# The natural variational form of the problem lives in $H^2(\Omega)$, +# +# ```math +# a(w,v) = \int_\Omega D^2 w : D^2 v, \qquad \ell(v) = \int_\Omega f\, v, +# ``` +# +# and an $H^2$-conforming space must be $C^1$ across element boundaries. That is +# expensive: the cheapest $C^1$ triangle is the quintic Argyris element, with 21 +# degrees of freedom (see its own tutorial). +# +# The **Morley element** buys simplicity by giving up conformity. It is the +# quadratic space $P_2(K)$ with only six degrees of freedom per triangle, +# +# ```math +# \ell^v(w) = w(v) \ \text{ at each vertex}, \qquad +# \ell^e(w) = \int_e \nabla w\cdot n \ \text{ on each edge}, +# ``` +# +# so neither $w$ nor $\nabla w$ is continuous across an edge — only the vertex +# values and the *mean* normal derivative are shared. The bilinear form is then +# applied element by element (a *broken* Hessian), and the method still converges. +# This is the classical nonconforming plate element of Morley. +# +# ## Discrete model +# +# We read an unstructured triangular mesh of the unit disk. The mesh is generated +# by `assets/nonstandard_reffes/MeshGenerator.jl` and carries a physical group +# named `"boundary"` on the rim. + +using Gridap +using GridapGmsh + +model = GmshDiscreteModel("../models/clamped_disk.msh") + +# ## FE space +# +# `morley` is a scalar element of order 2 on triangles. Its conformity is +# `H1Conformity()`, which in Gridap only says *which faces own and share degrees +# of freedom* — here vertices and edges. It is not a claim that the space is +# $C^0$; Morley is neither $C^0$ nor $C^1$. +# +# Both Morley degrees of freedom on the boundary — the vertex value and the mean +# normal derivative — are exactly the two clamped conditions, so the boundary +# condition is homogeneous. + +V = FESpace(model, ReferenceFE(morley, Float64, 2); dirichlet_tags="boundary") +U = TrialFESpace(V) + +# ## Weak form and solve +# +# `∇∇(u)` is the element-wise Hessian. On a conforming space it would be the +# distributional one; here the two differ, which is exactly the nonconformity. + +Ω = Triangulation(model) +dΩ = Measure(Ω, 6) + +f(x) = 1.0 + +a(u, v) = ∫( ∇∇(u) ⊙ ∇∇(v) )dΩ +l(v) = ∫( f * v )dΩ + +op = AffineFEOperator(a, l, U, V) +wh = solve(op) + +# ## Checking the result +# +# The exact deflection is $(1-r^2)^2/64$. We compare in the $L^2$ norm. The +# domain is meshed by straight-edged triangles, so the geometry itself is only +# approximated and the error will not go below that; on this mesh it is around a +# percent of the peak deflection. + +wex(x) = (1 - (x[1]^2 + x[2]^2))^2 / 64 + +wexh = CellField(wex, Ω) +e = wh - wexh +l2 = sqrt(sum( ∫( e * e )dΩ )) +peak = sqrt(sum( ∫( wexh * wexh )dΩ )) +println("relative L2 error: ", l2 / peak) + +# The deflection at the centre should be close to $1/64 = 0.015625$: + +println("exact centre deflection: ", 1 / 64) + +# ## Visualisation +# +# We write the deflection and the broken Hessian, whose jump across edges is what +# distinguishes a nonconforming method from a conforming one. + +mkpath("output_path") +writevtk(Ω, "output_path/morley", cellfields=["w" => wh, "hessian" => ∇∇(wh)]) + +# ## Nonconformity, made visible +# +# The value and the normal derivative do not match across an edge, but the *mean* +# of the normal derivative does — that is the shared degree of freedom. Averaged +# over each interior edge, the jump is at round-off: + +Λ = SkeletonTriangulation(model) +dΛ = Measure(Λ, 6) +n = get_normal_vector(Λ) + +facet_means = get_array( ∫( jump(∇(wh) ⋅ n) )dΛ ) +println("largest facet-mean jump of ∂w/∂n: ", maximum(abs, facet_means)) diff --git a/src/mtw_darcy_stokes.jl b/src/mtw_darcy_stokes.jl new file mode 100644 index 0000000..9b0dad8 --- /dev/null +++ b/src/mtw_darcy_stokes.jl @@ -0,0 +1,126 @@ +# In this tutorial, we will learn +# - How to discretise a problem that interpolates between Darcy and Stokes flow +# - How to use the Mardal--Tai--Winther element, which is stable uniformly in +# the parameter +# - Why an element can be $H(\operatorname{div})$-conforming and only *weakly* +# $H^1$ +# +# ## Problem statement +# +# The Darcy--Stokes (or Brinkman) problem: find the velocity $u$ and the pressure +# $p$ with +# +# ```math +# u - \varepsilon^2 \Delta u + \nabla p = f, \qquad \operatorname{div} u = 0, +# ``` +# +# on $\Omega = (0,1)^2$. At $\varepsilon \sim 1$ this is Stokes flow; as +# $\varepsilon \to 0$ it degenerates to Darcy flow. A discretisation that is +# stable for *both* limits, and uniformly in between, is the hard part: a Stokes +# pair typically degrades as $\varepsilon \to 0$, and a Darcy pair +# ($H(\operatorname{div})$-conforming, discontinuous tangentially) cannot control +# the viscous term. +# +# ## The element +# +# The **Mardal--Tai--Winther** element resolves this. On a triangle it is the +# 9-dimensional space +# +# ```math +# MTW(K) = P_1(K;\mathbb{R}^2) + \operatorname{curl}(b\,P_1(K)), +# ``` +# +# $b$ the cubic bubble, with three degrees of freedom per edge: two moments of +# $u\cdot n$ against $P_1(e)$, and one of $u\cdot t$ against the constant. +# +# Sharing the normal moments makes the space $H(\operatorname{div})$-conforming, +# so `div u = 0` is meaningful and the Darcy limit is fine. Sharing the single +# tangential moment makes it *weakly* $H^1$: the tangential trace jumps, but its +# mean over each edge does not, and that is enough to control the viscous term +# uniformly in $\varepsilon$. +# +# The same element exists on tetrahedra, where it is the element of Tai and +# Winther, with six degrees of freedom per face; `mtw` covers both. + +using Gridap +using Gridap.MultiField + +# ## Discrete model +# +# MTW is defined on simplices. + +n = 16 +model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) + +Ω = Triangulation(model) +dΩ = Measure(Ω, 8) + +# ## Manufactured solution +# +# The velocity is a curl, so it is divergence free by construction, and the +# stream function vanishes to second order on $\partial\Omega$, so $u$ vanishes +# there too. + +ψ(x) = (x[1]^2 * (1 - x[1])^2) * (x[2]^2 * (1 - x[2])^2) +uex(x) = VectorValue(∇(ψ)(x)[2], -∇(ψ)(x)[1]) +pex(x) = cos(pi * x[1]) * cos(pi * x[2]) + +# ## The solver +# +# The pressure is piecewise constant, with the mean fixed. Note `∇(u)` here is +# the *broken* gradient: the space is not $H^1$, and the viscous term is +# assembled element by element. + +function solve_darcy_stokes(ε) + f(x) = uex(x) - ε^2 * Δ(uex)(x) + ∇(pex)(x) + + V = FESpace(model, ReferenceFE(mtw, Float64, 1); dirichlet_tags="boundary") + Q = FESpace(model, ReferenceFE(lagrangian, Float64, 0); + conformity=:L2, constraint=:zeromean) + X = MultiFieldFESpace([TrialFESpace(V, uex), TrialFESpace(Q)]) + Y = MultiFieldFESpace([V, Q]) + + a((u, p), (v, q)) = ∫( u ⋅ v + ε^2 * (∇(u) ⊙ ∇(v)) - p * (∇ ⋅ v) - q * (∇ ⋅ u) )dΩ + l((v, q)) = ∫( f ⋅ v )dΩ + + solve(AffineFEOperator(a, l, X, Y)) +end + +# ## Uniform robustness +# +# The point of the element: the velocity error barely moves as $\varepsilon$ +# sweeps four orders of magnitude, from the Stokes regime to the Darcy one. + +for ε in (1.0, 1e-1, 1e-2, 1e-4) + uh, ph = solve_darcy_stokes(ε) + eu = uh - uex + ep = ph - pex + println("eps = ", ε, + " |u| = ", sqrt(sum( ∫( eu ⋅ eu )dΩ )), + " |p| = ", sqrt(sum( ∫( ep * ep )dΩ ))) +end + +# ## H(div)-conforming, weakly $H^1$ +# +# Interpolating a field outside the space shows both halves of the statement. The +# normal component is continuous outright; the tangential component jumps, but +# its mean over each edge is zero. + +w(x) = VectorValue(sin(3*x[1]) * x[2]^2, cos(2*x[2]) + x[1]^3) +V = FESpace(model, ReferenceFE(mtw, Float64, 1)) +wh = interpolate(w, V) + +Λ = SkeletonTriangulation(model) +dΛ = Measure(Λ, 8) +nΛ = get_normal_vector(Λ) +tΛ = Operation(n -> VectorValue(n[2], -n[1]))(nΛ) + +println("normal jump: ", sqrt(sum( ∫( jump(wh ⋅ nΛ) * jump(wh ⋅ nΛ) )dΛ ))) +println("tangential jump: ", sqrt(sum( ∫( jump(wh ⋅ tΛ) * jump(wh ⋅ tΛ) )dΛ ))) +println("mean tangential jump: ", maximum(abs, get_array( ∫( jump(wh ⋅ tΛ) )dΛ ))) + +# ## Visualisation + +uh, ph = solve_darcy_stokes(1e-4) +mkpath("output_path") +writevtk(Ω, "output_path/mtw", cellfields=["u" => uh, "p" => ph]) diff --git a/src/regge_metric.jl b/src/regge_metric.jl new file mode 100644 index 0000000..11ce6c5 --- /dev/null +++ b/src/regge_metric.jl @@ -0,0 +1,122 @@ +# In this tutorial, we will learn +# - How to use the Regge element, a symmetric-tensor-valued element that is +# *tangential--tangential* continuous +# - What its degrees of freedom mean geometrically +# - How it mirrors the Hellan--Herrmann--Johnson element +# +# ## The element +# +# The **Regge** element is the full symmetric-matrix-valued polynomial space +# $P_r(K;\mathbb{S})$, of any degree $r \ge 0$, with degrees of freedom +# +# ```math +# \int_e (t\cdot M t)\,\mu_i \ \text{ on each edge}, \qquad +# \int_K M : \tau \ \text{ in the cell}, +# ``` +# +# $t$ the unit tangent of the edge. Only the **tangential--tangential** component +# is shared between two triangles; everything else jumps. +# +# It is the exact mirror of the Hellan--Herrmann--Johnson element, which shares +# $n\cdot\sigma n$ instead — same space, same interior moments, same dimension, +# the normal swapped for the tangent. Their push-forwards differ accordingly +# (double covariant rather than double contravariant Piola), and both end up with +# the same, diagonal, change of basis. +# +# ## What the degrees of freedom mean +# +# Regge elements come from Regge calculus, where a geometry is described not by +# coordinates but by **edge lengths**. That is exactly what the lowest-order +# degrees of freedom are. If $M$ is a perturbation of the metric, the induced +# change in the squared length of an edge $e$ with tangent $t$ and length $|e|$ is +# +# ```math +# \delta(|e|^2) = |e| \int_e t\cdot M t , +# ``` +# +# i.e. the Regge degree of freedom, up to the edge length. A discrete metric in +# the lowest-order Regge space is *precisely* an assignment of one number per +# edge — which is why the tangential--tangential component, and only that, has to +# be single-valued: two triangles sharing an edge must agree on its length. + +using Gridap +using Gridap.TensorValues +using Gridap.Geometry + +# ## Discrete model + +n = 8 +model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) + +Ω = Triangulation(model) +dΩ = Measure(Ω, 6) + +# ## The space +# +# `regge` is a symmetric-tensor-valued element of degree `r`. At `r = 0` it has +# exactly one degree of freedom per edge and none in the cell. + +r = 0 +V = FESpace(model, ReferenceFE(regge, Float64, r)) + +topo = get_grid_topology(model) +println("dofs: ", num_free_dofs(V), " edges: ", num_faces(topo, 1)) + +# ## Interpolation +# +# The space is all of $P_r(K;\mathbb{S})$, so it reproduces any symmetric tensor +# field of that degree exactly. Here `r = 0`: constants. + +Mconst(x) = SymTensorValue{2,Float64}(2.0, -1.0, 3.0) +Mh = interpolate(Mconst, V) +e = Mh - Mconst +println("interpolation error, constant metric: ", sqrt(sum( ∫( e ⊙ e )dΩ ))) + +# ## Tangential--tangential continuity +# +# The defining property. We interpolate a field that is *not* in the space, so +# the jumps are not zero for a trivial reason, and measure the jump of +# $t\cdot M t$ against the jump of the whole tensor. + +g(x) = SymTensorValue{2,Float64}(sin(2*x[1]), cos(3*x[2]), sin(x[1] + x[2])) +gh = interpolate(g, V) + +Λ = SkeletonTriangulation(model) +dΛ = Measure(Λ, 6) +nΛ = get_normal_vector(Λ).⁺ +tΛ = Operation(n -> VectorValue(n[2], -n[1]))(nΛ) + +tt = sqrt(sum( ∫( (tΛ ⋅ jump(gh) ⋅ tΛ) * (tΛ ⋅ jump(gh) ⋅ tΛ) )dΛ )) +full = sqrt(sum( ∫( jump(gh) ⊙ jump(gh) )dΛ )) +println("jump of the tt component: ", tt) +println("jump of the full tensor: ", full) + +# Compare with the mirror element, which shares the normal--normal component +# instead — on the same mesh, with the same field, the two continuities are +# exchanged: + +Vh = FESpace(model, ReferenceFE(hhj, Float64, r)) +hh = interpolate(g, Vh) +println("HHJ, jump of the nn component: ", + sqrt(sum( ∫( (nΛ ⋅ jump(hh) ⋅ nΛ) * (nΛ ⋅ jump(hh) ⋅ nΛ) )dΛ ))) +println("HHJ, jump of the tt component: ", + sqrt(sum( ∫( (tΛ ⋅ jump(hh) ⋅ tΛ) * (tΛ ⋅ jump(hh) ⋅ tΛ) )dΛ ))) + +# ## Higher degree +# +# The element is a family. At degree `r` there are `r+1` degrees of freedom per +# edge and `3r(r+1)/2` in the cell, and the space reproduces any +# $P_r(K;\mathbb{S})$ field exactly. + +for r in 0:2 + Vr = FESpace(model, ReferenceFE(regge, Float64, r)) + Mr(x) = SymTensorValue{2,Float64}(1.0 + x[1]^r, 2.0 - x[2]^r, 3.0 + (x[1]*x[2])^(r ÷ 2)) + er = interpolate(Mr, Vr) - Mr + println("r = ", r, " dofs = ", num_free_dofs(Vr), + " interpolation error = ", sqrt(sum( ∫( er ⊙ er )dΩ ))) +end + +# ## Visualisation + +mkpath("output_path") +writevtk(Ω, "output_path/regge", cellfields=["M" => gh]) From 6205bf8cdadcea3870a982d470f32129c9073b17 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Thu, 10 Sep 2026 09:08:19 +1000 Subject: [PATCH 2/4] Activate tutorials --- deps/build.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/deps/build.jl b/deps/build.jl index 1f7f8c5..3cb7adb 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -30,6 +30,13 @@ files = [ "Poisson with HHO on polytopal meshes"=>"poisson_hho.jl", "Block assembly and solvers: Incompressible Stokes example"=>"stokes_blocks.jl", "Lagrange multipliers" => "lagrange_multipliers.jl", + "Biharmonic equation (with Morley)" => "morley_biharmonic.jl", + "Biharmonic equation (with Argyris)" => "argyris_biharmonic.jl", + "Kirchhoff plate (with HHJ)" => "hhj_plate.jl", + "Mixed elasticity (with Arnold-Winther)" => "arnold_winther_elasticity.jl", + "Darcy-Stokes (with Mardal-Tai-Winther)" => "mtw_darcy_stokes.jl", + "Regge metrics" => "regge_metric.jl", + "Deviatoric stress (with GLS)" => "gls_stokes.jl", "Low-level API - Arrays"=>"arrays_dev.jl", "Low-level API - Geometry" => "geometry_dev.jl", "Low-level API - ReferenceFEs" => "reffes_dev.jl", From fd8e589ae5ae386042c1249dff7a897acddfe6b3 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 16 Sep 2026 17:16:05 +1000 Subject: [PATCH 3/4] Added more tutorials --- deps/build.jl | 3 + src/hermite_beam.jl | 282 ++++++++++++++++++++++++++++++++ src/laminate_plate.jl | 229 ++++++++++++++++++++++++++ src/weak_symmetry_elasticity.jl | 185 +++++++++++++++++++++ 4 files changed, 699 insertions(+) create mode 100644 src/hermite_beam.jl create mode 100644 src/laminate_plate.jl create mode 100644 src/weak_symmetry_elasticity.jl diff --git a/deps/build.jl b/deps/build.jl index 3cb7adb..720715c 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -30,10 +30,13 @@ files = [ "Poisson with HHO on polytopal meshes"=>"poisson_hho.jl", "Block assembly and solvers: Incompressible Stokes example"=>"stokes_blocks.jl", "Lagrange multipliers" => "lagrange_multipliers.jl", + "Euler-Bernoulli beam (with Hermite)" => "hermite_beam.jl", "Biharmonic equation (with Morley)" => "morley_biharmonic.jl", "Biharmonic equation (with Argyris)" => "argyris_biharmonic.jl", "Kirchhoff plate (with HHJ)" => "hhj_plate.jl", "Mixed elasticity (with Arnold-Winther)" => "arnold_winther_elasticity.jl", + "Elasticity with weak symmetry (with stacked BDM)" => "weak_symmetry_elasticity.jl", + "Laminate plate (with a Lagrange-Argyris product)" => "laminate_plate.jl", "Darcy-Stokes (with Mardal-Tai-Winther)" => "mtw_darcy_stokes.jl", "Regge metrics" => "regge_metric.jl", "Deviatoric stress (with GLS)" => "gls_stokes.jl", diff --git a/src/hermite_beam.jl b/src/hermite_beam.jl new file mode 100644 index 0000000..c14d89a --- /dev/null +++ b/src/hermite_beam.jl @@ -0,0 +1,282 @@ +# In this tutorial, we will learn +# - How to use the Hermite element, the cubic simplex that carries vertex +# gradients as degrees of freedom +# - How to solve the Euler--Bernoulli beam, the one setting where that element +# is $C^1$ and the fourth-order problem is therefore conforming +# - Why the same element is only $C^0$ in two dimensions, and what it is worth +# there +# +# ## The element +# +# The **Hermite** element [Ciarlet & Raviart, *Arch. Rational Mech. Anal.* 46 +# (1972) 177] is the full cubic space $P_3(K)$ on a simplex +# $K \subset \mathbb{R}^D$, with degrees of freedom +# +# ```math +# \ell^{v,0}(u) = u(v), \qquad +# \ell^{v,i}(u) = \partial_i u(v), \quad i = 1,\dots,D, +# ``` +# +# at each vertex $v$, plus the value $u(b)$ at the barycentre $b$ of each +# two-dimensional face. That is $D+1$ degrees of freedom per vertex and one per +# triangular face: 4 on a segment, 10 on a triangle, 20 on a tetrahedron — in +# each case exactly $\dim P_3(K)$. +# +# Because a vertex owns its *whole* gradient, and that gradient is stated in the +# global Cartesian frame, every cell meeting at a vertex shares it. The global +# space is therefore $C^0$ **and** has a single-valued gradient at every vertex. +# It is *not* $C^1$: across an edge the gradients of the two neighbouring cells +# agree at the two endpoints and nowhere else. +# +# Like Morley and Argyris, Hermite is mapped by the plain pullback, and its +# gradient degrees of freedom are not preserved by it — $\nabla u(v) = J^{-T} +# \hat\nabla\hat u(\hat v)$ — so a per-cell change of basis is applied. That is +# internal: `ReferenceFE(hermite, Float64)` is all a user writes. The element +# exists for order 3 only, so the order argument is optional. + +using Gridap +using Gridap.ReferenceFEs + +# ## The Euler--Bernoulli beam +# +# In one dimension the vertices are the *only* interfaces, so "$C^1$ at the +# vertices" is $C^1$ outright and the element is $H^2$-conforming. This is the +# classical beam element of every structural analysis textbook. +# +# We solve the **clamped beam** under a uniform load: find the deflection $w$ of +# a beam of unit length and unit bending stiffness with +# +# ```math +# w'''' = f \ \text{ in } (0,1), \qquad +# w = w' = 0 \ \text{ at } x = 0, 1 , +# ``` +# +# whose variational form is +# +# ```math +# a(w,v) = \int_0^1 w'' v'' , \qquad \ell(v) = \int_0^1 f\, v . +# ``` +# +# For $f \equiv 1$ the deflection is the quartic +# +# ```math +# w(x) = \frac{x^2 (1-x)^2}{24}, +# ``` +# +# so the midspan deflection is exactly $1/384$. + +f = 1.0 +w(x) = x[1]^2 * (1 - x[1])^2 / 24 + +# ## Discrete model and space +# +# A one-dimensional Cartesian model has `SEGMENT` cells, and its face labeling +# already provides a `"boundary"` tag for the two endpoints. +# +# The clamped condition is where this element pays off. A boundary vertex owns +# exactly two degrees of freedom, $w(v)$ and $w'(v)$, and those are precisely +# the two clamped conditions — so `dirichlet_tags="boundary"` imposes the whole +# boundary condition, with nothing left over and nothing missing. The data is +# homogeneous, so `TrialFESpace(V)` suffices. + +model = CartesianDiscreteModel((0, 1), (8,)) + +V = FESpace(model, ReferenceFE(hermite, Float64); dirichlet_tags="boundary") +U = TrialFESpace(V) + +reffe = ReferenceFE(SEGMENT, hermite, Float64) +println("dofs per cell: ", num_dofs(reffe)) +println("free dofs: ", num_free_dofs(V), " clamped: ", num_dirichlet_dofs(V)) + +# ## Weak form and solve +# +# `∇∇(u)` is the Hessian, a $1\times1$ tensor in one dimension, so the bilinear +# form is the double contraction. Unlike the Morley and HHJ tutorials this is +# *not* a broken form: the space really is $H^2$, so the element-wise Hessian +# and the distributional one agree. + +Ω = Triangulation(model) +dΩ = Measure(Ω, 8) + +a(u, v) = ∫( ∇∇(u) ⊙ ∇∇(v) )dΩ +l(v) = ∫( f * v )dΩ + +wh = solve(AffineFEOperator(a, l, U, V)) + +e = w - wh +println("L2 error: ", sqrt(sum( ∫( e * e )dΩ ))) +println("exact midspan deflection: ", 1 / 384) + +# ## Convergence +# +# The exact deflection is quartic, so it is not reproduced exactly and we see +# genuine rates: $O(h^4)$ in $L^2$ and $O(h^2)$ in the $H^2$ seminorm, which is +# optimal for a cubic conforming space. + +function beam(n) + model = CartesianDiscreteModel((0, 1), (n,)) + V = FESpace(model, ReferenceFE(hermite, Float64); dirichlet_tags="boundary") + U = TrialFESpace(V) + Ω = Triangulation(model) + dΩ = Measure(Ω, 8) + a(u, v) = ∫( ∇∇(u) ⊙ ∇∇(v) )dΩ + l(v) = ∫( f * v )dΩ + wh = solve(AffineFEOperator(a, l, U, V)) + e = w - wh + el2 = sqrt(sum( ∫( e * e )dΩ )) + eh2 = sqrt(sum( ∫( ∇∇(e) ⊙ ∇∇(e) )dΩ )) + el2, eh2 +end + +for n in (2, 4, 8, 16) + el2, eh2 = beam(n) + println("n = ", n, " L2 = ", el2, " H2 = ", eh2) +end + +# ## The same element in two dimensions +# +# On triangles the edges are interfaces too, and the element is only $C^0$. The +# fourth-order problem is out of reach — that is what the Morley and Argyris +# tutorials are for — but the vertex gradients are still shared, and that makes +# the global space a proper *subspace* of the continuous cubic Lagrange space. +# +# The two count differently. On a triangulation with $V$ vertices, $E$ edges and +# $C$ cells, cubic Lagrange has one degree of freedom per vertex, two per edge +# and one per cell, while Hermite has three per vertex and one per cell: +# +# ```math +# \dim P_3^{\text{Lagrange}} = V + 2E + C, \qquad +# \dim P_3^{\text{Hermite}} = 3V + C . +# ``` +# +# On a large triangulation $E \approx 3V$ and $C \approx 2V$, so Hermite is +# about $5V$ against about $9V$ — a little over half the degrees of freedom for +# the same local polynomial space. + +using Gridap.Geometry + +model2 = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (8, 8))) +topo = get_grid_topology(model2) +nV, nE, nC = num_faces(topo, 0), num_faces(topo, 1), num_faces(topo, 2) + +Vh = FESpace(model2, ReferenceFE(hermite, Float64)) +Vl = FESpace(model2, ReferenceFE(lagrangian, Float64, 3)) + +println("hermite: ", num_free_dofs(Vh), " 3V + C = ", 3*nV + nC) +println("lagrange: ", num_free_dofs(Vl), " V + 2E + C = ", nV + 2*nE + nC) + +# ## A cheaper cubic +# +# We solve a reaction--diffusion problem with natural boundary conditions, +# +# ```math +# -\Delta u + u = f \ \text{ in } \Omega = (0,1)^2, \qquad +# \frac{\partial u}{\partial n} = 0 \ \text{ on } \partial\Omega, +# ``` +# +# and compare the two spaces on the same meshes. Natural conditions keep the +# comparison clean: a Hermite vertex on the boundary owns its full gradient, so +# `dirichlet_tags` there would pin the *normal* derivative as well as the value, +# which plain Dirichlet data does not determine. (In the beam that was a +# feature — clamping *is* a condition on $w$ and $w'$.) +# +# Both spaces converge at $O(h^4)$ in $L^2$ and $O(h^3)$ in $H^1$ — on these +# coarse meshes the Hermite rates are still climbing towards those values, the +# Lagrange ones are there already. Hermite carries a somewhat larger error +# constant, and about 57% of the degrees of freedom. + +uex(x) = cos(pi*x[1]) * cos(pi*x[2]) + cos(2*pi*x[1]) +fex(x) = -Δ(uex)(x) + uex(x) + +function reaction_diffusion(n, reffe) + model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) + V = FESpace(model, reffe) + Ω = Triangulation(model) + dΩ = Measure(Ω, 8) + a(u, v) = ∫( ∇(u) ⋅ ∇(v) + u * v )dΩ + l(v) = ∫( fex * v )dΩ + uh = solve(AffineFEOperator(a, l, V, V)) + e = uex - uh + el2 = sqrt(sum( ∫( e * e )dΩ )) + eh1 = sqrt(sum( ∫( ∇(e) ⋅ ∇(e) )dΩ )) + num_free_dofs(V), el2, eh1 +end + +for n in (2, 4, 8, 16) + ndh, l2h, h1h = reaction_diffusion(n, ReferenceFE(hermite, Float64)) + ndl, l2l, h1l = reaction_diffusion(n, ReferenceFE(lagrangian, Float64, 3)) + println("n = ", n, + " hermite: ", ndh, " dofs, L2 = ", l2h, ", H1 = ", h1h, + " lagrange: ", ndl, " dofs, L2 = ", l2l, ", H1 = ", h1l) +end + +# ## $C^0$, and $C^1$ at the vertices only +# +# The claim about the space is worth checking rather than asserting. We +# interpolate a smooth function and look at the interpolant from inside each +# cell separately. +# +# Evaluating a `CellField` at a `CellPoint` built from the reference vertices +# gives, for every cell, the gradient at its own three vertices. Gathering those +# per mesh vertex, the *spread* over the cells meeting there is zero for Hermite +# and $O(h^3)$ for Lagrange: the Hermite gradient at a vertex is a single shared +# degree of freedom, the Lagrange one is whatever each cell happens to produce. + +using Gridap.CellData +using FillArrays +using LinearAlgebra + +g(x) = cos(pi*x[1]) * cos(pi*x[2]) + cos(2*pi*x[1]) + +Ω2 = Triangulation(model2) +cell_to_vertices = get_faces(topo, 2, 0) +cell_points = CellPoint( + Fill(get_vertex_coordinates(TRI), num_cells(model2)), Ω2, ReferenceDomain() +) + +Λ = SkeletonTriangulation(model2) +dΛ = Measure(Λ, 8) + +for (name, space) in (("hermite ", Vh), ("lagrange", Vl)) + gh = interpolate(g, space) + cell_vertex_grads = ∇(gh)(cell_points) + + at_vertex = [VectorValue{2,Float64}[] for _ in 1:nV] + for cell in 1:num_cells(model2) + for (i, vertex) in enumerate(cell_to_vertices[cell]) + push!(at_vertex[vertex], cell_vertex_grads[cell][i]) + end + end + spread = maximum(at_vertex) do vals + maximum(v -> norm(v - first(vals)), vals) + end + + jump_u = sqrt(sum( ∫( jump(gh) * jump(gh) )dΛ )) + jump_g = sqrt(sum( ∫( jump(∇(gh)) ⋅ jump(∇(gh)) )dΛ )) + println(name, " vertex gradient spread: ", spread, + " [[u]]: ", jump_u, " [[grad u]]: ", jump_g) +end + +# The skeleton jumps say the rest: the value jump is at round-off for both +# spaces — they are $C^0$ — while the gradient jump is $O(1)$ for both. Hermite +# is not $C^1$, and its gradient jump is in fact the larger of the two here. The +# vertex gradients are shared, not the edge ones. +# +# For the Hermite interpolant that shared value is also the *exact* one, since +# the degree of freedom is $\partial_i u(v)$ itself: + +gh = interpolate(g, Vh) +cell_vertex_grads = ∇(gh)(cell_points) +vertex_coords = get_vertex_coordinates(topo) + +err = maximum( + norm(cell_vertex_grads[cell][i] - ∇(g)(vertex_coords[vertex])) + for cell in 1:num_cells(model2) + for (i, vertex) in enumerate(cell_to_vertices[cell]) +) +println("largest vertex gradient error of the interpolant: ", err) + +# ## Visualisation + +mkpath("output_path") +writevtk(Ω2, "output_path/hermite", cellfields=["u" => gh, "grad" => ∇(gh)]) diff --git a/src/laminate_plate.jl b/src/laminate_plate.jl new file mode 100644 index 0000000..7aa399f --- /dev/null +++ b/src/laminate_plate.jl @@ -0,0 +1,229 @@ +# In this tutorial, we will learn +# - How to build a single element whose components have *different* +# conformities, with `CartProdRefFE` +# - How to solve a plate problem in which stretching and bending are coupled, +# using one displacement field rather than two unknowns +# - How to write kinematics as slices of $\nabla U$ and $\nabla\nabla U$ +# +# ## Problem statement +# +# A **laminated plate** is a stack of bonded layers. If the stack is symmetric +# about the mid-surface, stretching and bending do not interact and the two +# problems are solved independently — that is the usual situation, and it is why +# plate tutorials normally treat bending alone. If the stack is *unsymmetric*, +# with stiffer layers above than below, the two are coupled: pulling on it +# bends it, and bending it stretches it. +# +# Classical laminate theory writes this with three constitutive tensors, $A$ for +# stretching, $D$ for bending and $B$ for the coupling, $B$ vanishing exactly +# when the stack is symmetric. We take the simplest instance of it: all three +# proportional to the identity, with zero Poisson ratio, so that a single number +# $b$ carries the coupling. The unknowns are the in-plane displacement +# $u = (u_1,u_2)$ and the transverse deflection $w$, and the energy is +# +# ```math +# E(u,w) = \frac{1}{2}\int_\Omega +# \varepsilon(u) : \varepsilon(u) +# + 2b\, \varepsilon(u) : \kappa(w) +# + \kappa(w) : \kappa(w) +# \;-\; \int_\Omega f \cdot (u_1,u_2,w) , +# ``` +# +# with $\varepsilon(u)$ the symmetric gradient and $\kappa(w) = \nabla\nabla w$ +# the curvature. The energy is positive definite for $|b| < 1$. At $b = 0$ it is +# the sum of a linear elasticity problem and a Kirchhoff plate problem. +# +# The two terms ask for different smoothness: $\varepsilon(u)$ needs +# $u \in H^1$, while $\kappa(w)$ needs $w \in H^2$. So the discrete space must +# be $C^0$ in its first two components and $C^1$ in the third. +# +# ## One field, not two +# +# Nothing *forces* the three components into one element — a `MultiFieldFESpace` +# pairing a vector Lagrange space with an Argyris space describes the same +# discrete problem. But these are not three unknowns that happen to be solved +# together; they are the three components of one displacement. If you +# never need them apart, the multi-field machinery is bookkeeping you are paying +# for and not using: two spaces to build, blocked trial and test functions to +# unpack in every term, a solution that comes back in pieces. +# +# `CartProdRefFE(r₁, r₂, r₃)` builds the product element instead. The factors +# must share a polytope and a value type — that is what makes the result one +# `MultiValue` rather than a multi-field space — but they need share nothing +# else, neither degree nor conformity. Here all three factors are scalar +# elements on a triangle, so the product is `VectorValue{3}`-valued, with +# factor $c$ occupying component $c$ of the value and block $c$ of the degrees +# of freedom. + +using Gridap +using Gridap.ReferenceFEs +using Gridap.TensorValues + +# ## Discrete model and space +# +# Argyris is defined on triangles, so we simplexify a Cartesian mesh. +# +# The membrane components are quadratic Lagrange and the bending component is +# quintic Argyris. Their conformities differ — `H1` against `H2` — so the +# product reports a `CartProdConformity`, which applies each factor's own +# ownership rules to its own block. That is the whole mechanism: the first two +# components glue like Lagrange, the third like Argyris. + +n = 8 +model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) + +membrane_reffe = ReferenceFE(TRI, lagrangian, Float64, 2) +bending_reffe = ReferenceFE(TRI, argyris, Float64, 5) +reffe = CartProdRefFE(membrane_reffe, membrane_reffe, bending_reffe) + +println("dofs per cell: ", num_dofs(reffe), + " = ", num_dofs(membrane_reffe), " + ", num_dofs(membrane_reffe), + " + ", num_dofs(bending_reffe)) +println("conformity: ", Conformity(reffe)) + +# The plate is clamped all round: $u = 0$ and $w = \partial_n w = 0$ on +# $\partial\Omega$. Setting every boundary degree of freedom to zero also pins +# the vertex Hessians of $w$, which is slightly more than a clamp requires — the +# Argyris tutorial avoids this by prescribing a manufactured solution on the +# boundary. Here we keep the homogeneous data, and note that every comparison +# below is between discretisations that share it. + +V = FESpace(model, reffe; dirichlet_tags="boundary") +U = TrialFESpace(V) + +Ω = Triangulation(model) +dΩ = Measure(Ω, 12) + +# ## Kinematics +# +# This is where one field pays. The unknown `U` is `VectorValue{3}`-valued on a +# two-dimensional domain, so `∇(U)` is $2\times3$ with +# $(\nabla U)_{ij} = \partial_i U_j$, and `∇∇(U)` is $2\times2\times3$. Both the +# membrane strain and the curvature are slices of those two objects, obtained by +# contracting the *component* index: +# +# ```math +# \varepsilon(u) = \operatorname{sym}\big( \nabla U \cdot P \big), \qquad +# \kappa(w) = \nabla\nabla U \cdot e_3 , +# ``` +# +# where $P$ keeps the first two components and $e_3$ the third. Contraction with +# `⋅` acts on the last index of its left argument, which is the component index +# in both cases, so the derivative indices are untouched. + +P = TensorValue{3,2}(1.0, 0.0, 0.0, 0.0, 1.0, 0.0) +e3 = VectorValue(0.0, 0.0, 1.0) + +displacement(U) = U ⋅ P +deflection(U) = U ⋅ e3 +strain(U) = symmetric_part( ∇(U) ⋅ P ) +curvature(U) = ∇∇(U) ⋅ e3 + +# ## Weak form and solve +# +# The load is a uniform body force with an in-plane and a transverse part. It +# enters as `∫( f ⋅ v )`, a single contraction against the whole test function — +# there is no unpacking to do. + +f = VectorValue(1.0, 1.0, 1.0) + +function solve_laminate(b) + a(u, v) = ∫( strain(u) ⊙ strain(v) + + b * ( strain(u) ⊙ curvature(v) + curvature(u) ⊙ strain(v) ) + + curvature(u) ⊙ curvature(v) )dΩ + l(v) = ∫( f ⋅ v )dΩ + solve(AffineFEOperator(a, l, U, V)) +end + +Uh = solve_laminate(0.5) + +println("compliance: ", sum( ∫( f ⋅ Uh )dΩ )) + +# ## At $b = 0$ the product is exactly the two spaces +# +# With the coupling switched off the energy splits, and the product space must +# reproduce, component by component, what the two spaces give on their own. We +# solve the same two problems in a plain vector Lagrange space and a plain +# Argyris space and compare. The degrees of freedom add up exactly, and so do +# the solutions. + +U0 = solve_laminate(0.0) + +Vw = FESpace(model, ReferenceFE(argyris, Float64, 5); dirichlet_tags="boundary") +aw(w, v) = ∫( ∇∇(w) ⊙ ∇∇(v) )dΩ +lw(v) = ∫( 1.0 * v )dΩ +w_solo = solve(AffineFEOperator(aw, lw, TrialFESpace(Vw), Vw)) + +Vu = FESpace(model, ReferenceFE(lagrangian, VectorValue{2,Float64}, 2); + dirichlet_tags="boundary") +au(u, v) = ∫( ε(u) ⊙ ε(v) )dΩ +lu(v) = ∫( VectorValue(1.0, 1.0) ⋅ v )dΩ +u_solo = solve(AffineFEOperator(au, lu, TrialFESpace(Vu), Vu)) + +println("free dofs: ", num_free_dofs(V), + " = ", num_free_dofs(Vu), " + ", num_free_dofs(Vw)) + +ew = deflection(U0) - w_solo +eu = displacement(U0) - u_solo +println("deflection difference: ", sqrt(sum( ∫( ew * ew )dΩ ))) +println("displacement difference: ", sqrt(sum( ∫( eu ⋅ eu )dΩ ))) + +# ## The coupling +# +# Now load the plate transversely only, with $f = (0,0,1)$, and vary $b$. A +# symmetric laminate does not stretch under a transverse load, so the in-plane +# displacement is *exactly* zero at $b = 0$ — the two blocks of the stiffness +# matrix do not talk to each other, and the in-plane load is what was removed. +# As soon as $b \neq 0$ it is not, and the plate also deflects more: the +# stretching it is now free to do lowers the energy at fixed load. + +f_transverse = VectorValue(0.0, 0.0, 1.0) + +for b in (0.0, 0.25, 0.5, 0.75) + a(u, v) = ∫( strain(u) ⊙ strain(v) + + b * ( strain(u) ⊙ curvature(v) + curvature(u) ⊙ strain(v) ) + + curvature(u) ⊙ curvature(v) )dΩ + l(v) = ∫( f_transverse ⋅ v )dΩ + Ub = solve(AffineFEOperator(a, l, U, V)) + nu = sqrt(sum( ∫( displacement(Ub) ⋅ displacement(Ub) )dΩ )) + nw = sqrt(sum( ∫( deflection(Ub) * deflection(Ub) )dΩ )) + println("b = ", b, " |u| = ", nu, " |w| = ", nw, + " compliance = ", sum( ∫( f_transverse ⋅ Ub )dΩ )) +end + +# ## One field, two conformities +# +# The claim to check is that the components really do glue differently. We +# interpolate a smooth function into the *unconstrained* product space and look +# at the jumps across the skeleton, component by component. +# +# The whole field is $C^0$, as both factors are. But the gradient jumps in the +# membrane components, which are only Lagrange, and vanishes in the bending +# component, which is Argyris. One element, two continuities. + +Vfree = FESpace(model, reffe) +g(x) = VectorValue(sin(2*x[1]) * x[2], cos(3*x[2]), sin(2*x[1]) * cos(3*x[2])) +gh = interpolate(g, Vfree) + +Λ = SkeletonTriangulation(model) +dΛ = Measure(Λ, 12) + +jump_u = jump(∇(gh)) ⋅ P +jump_w = jump(∇(gh)) ⋅ e3 + +println("value jump, whole field: ", + sqrt(sum( ∫( jump(gh) ⋅ jump(gh) )dΛ ))) +println("gradient jump, membrane (P2): ", + sqrt(sum( ∫( jump_u ⊙ jump_u )dΛ ))) +println("gradient jump, bending (Arg): ", + sqrt(sum( ∫( jump_w ⋅ jump_w )dΛ ))) + +# ## Visualisation +# +# The solution is a single `FEFunction`, so it is written out as one field. + +mkpath("output_path") +writevtk(Ω, "output_path/laminate_plate", + cellfields=["U" => Uh, + "u" => displacement(Uh), + "w" => deflection(Uh)]) diff --git a/src/weak_symmetry_elasticity.jl b/src/weak_symmetry_elasticity.jl new file mode 100644 index 0000000..3ddb09a --- /dev/null +++ b/src/weak_symmetry_elasticity.jl @@ -0,0 +1,185 @@ +# In this tutorial, we will learn +# - How to solve linear elasticity in mixed form with *weakly* imposed symmetry +# - How to build a tensor-valued $H(\mathrm{div})$ space by stacking copies of a +# vector-valued element with `CartProdRefFE` +# - Why the copies are stacked on the *last* index, and what that changes in the +# weak form +# +# ## Problem statement +# +# The Arnold--Winther tutorial solves the Hellinger--Reissner system with +# a stress space that is symmetric pointwise. Building such a space is hard. The +# alternative is to drop the symmetry from the space and impose it with a +# Lagrange multiplier $\gamma$, the *rotation*. With the identity compliance the +# system is: find $(\sigma, u, \gamma)$ with +# +# ```math +# \begin{aligned} +# \int_\Omega \sigma : \tau + \int_\Omega u\cdot\operatorname{div}\tau +# + \int_\Omega \gamma\, (\chi : \tau) &= 0, \\ +# \int_\Omega v\cdot\operatorname{div}\sigma &= -\int_\Omega f\cdot v, \\ +# \int_\Omega \eta\, (\chi : \sigma) &= 0, +# \end{aligned} +# ``` +# +# for all $(\tau, v, \eta)$, where +# +# ```math +# \chi = \begin{pmatrix} 0 & 1 \\ -1 & 0\end{pmatrix}, +# \qquad \chi : \tau = \tau_{12} - \tau_{21}, +# ``` +# +# so the third equation says exactly that $\sigma$ is symmetric in the weak +# sense, and $\gamma$ converges to $\tfrac12(\partial_1 u_2 - \partial_2 u_1)$. +# The stress is now an arbitrary tensor field with rows -- or columns, see below +# -- in $H(\operatorname{div})$, which is a space we already have. +# +# ## The stress space +# +# `CartProdRefFE(reffe, K)` stacks `K` independent copies of any reference +# element. With a vector-valued atom the values become tensors, +# +# ```math +# \texttt{VectorValue\{d\}} \;\longmapsto\; \texttt{TensorValue\{d,K\}}, +# ``` +# +# one copy per **column**. Two copies of a $2$-dimensional $H(\operatorname{div})$ +# element therefore give a tensor field whose every column is +# $H(\operatorname{div})$ -- the stress space of elasticity with weakly imposed +# symmetry. +# +# The copies go on the *last* index because Gridap writes the derivative index +# *first*, $(\nabla A)_{kij} = \partial_k A_{ij}$. Appending the copy index is +# then the one placement that commutes with differentiation, and one consequence +# is exactly what this tutorial needs: since +# $\operatorname{div} A = \operatorname{tr}(\nabla A)$ and `tr` of a third-order +# tensor traces its first two indices, +# +# ```math +# (\operatorname{div} A)_j = \partial_i A_{ij}, +# ``` +# +# so `∇⋅S` on a stacked space is precisely the vector of the copies' +# divergences. The literature writes this element row-wise instead; the system +# above is its transpose, and since the exact stress is symmetric the two +# describe the same solution. +# +# ## The elements +# +# Arnold, Falk and Winther's family on simplices takes, for $k \ge 1$, +# +# ```math +# \Sigma_h = \mathrm{BDM}_k^{\,2}, \qquad +# V_h = P_{k-1}^{\mathrm{dc}}(\mathbb{R}^2), \qquad +# Q_h = P_{k-1}^{\mathrm{dc}}, +# ``` +# +# and converges at $O(h^k)$ in all three unknowns. `CartProdRefFE` builds +# $\Sigma_h$ from the `bdm` element already in Gridap. + +using Gridap +using Gridap.TensorValues +using Gridap.MultiField + +# ## Discrete model + +n = 8 +model = simplexify(CartesianDiscreteModel((0, 1, 0, 1), (n, n))) + +Ω = Triangulation(model) +degree = 8 +dΩ = Measure(Ω, degree) + +# ## Manufactured solution +# +# The same displacement as the Arnold--Winther tutorial, vanishing on +# $\partial\Omega$, together with its stress $\sigma = \varepsilon(u)$, its +# rotation $\gamma$, and the body force. Note $\sigma$ is written as a full +# `TensorValue`: the discrete stress is not symmetric, so the error cannot be +# measured against a `SymTensorValue`. + +const π2 = pi^2 +uex(x) = VectorValue(sin(pi*x[1]) * sin(pi*x[2]), sin(2*pi*x[1]) * sin(pi*x[2])) + +function σex(x) + u1x = pi * cos(pi*x[1]) * sin(pi*x[2]) + u1y = pi * sin(pi*x[1]) * cos(pi*x[2]) + u2x = 2*pi * cos(2*pi*x[1]) * sin(pi*x[2]) + u2y = pi * sin(2*pi*x[1]) * cos(pi*x[2]) + s12 = (u1y + u2x) / 2 + TensorValue{2,2,Float64}(u1x, s12, s12, u2y) +end + +function γex(x) + u1y = pi * sin(pi*x[1]) * cos(pi*x[2]) + u2x = 2*pi * cos(2*pi*x[1]) * sin(pi*x[2]) + (u2x - u1y) / 2 +end + +function fex(x) + u1 = sin(pi*x[1]) * sin(pi*x[2]) + u2 = sin(2*pi*x[1]) * sin(pi*x[2]) + d11u1, d22u1 = -π2*u1, -π2*u1 + d12u1 = π2 * cos(pi*x[1]) * cos(pi*x[2]) + d11u2, d22u2 = -4*π2*u2, -π2*u2 + d12u2 = 2*π2 * cos(2*pi*x[1]) * cos(pi*x[2]) + VectorValue(-(d11u1 + (d22u1 + d12u2)/2), -((d12u1 + d11u2)/2 + d22u2)) +end + +# ## The solver +# +# The three spaces, and the system written exactly as above. The displacement +# boundary condition $u = 0$ is natural here, so no `dirichlet_tags` appear, and +# no integration by parts is needed: `∇⋅S` is available directly on the stacked +# space. + +const χ = TensorValue{2,2,Float64}(0.0, -1.0, 1.0, 0.0) + +function solve_weak_symmetry(k) + Σ = FESpace(model, CartProdRefFE(ReferenceFE(TRI, bdm, Float64, k), 2)) + V = FESpace(model, ReferenceFE(lagrangian, VectorValue{2,Float64}, k-1); conformity=:L2) + Q = FESpace(model, ReferenceFE(lagrangian, Float64, k-1); conformity=:L2) + X = MultiFieldFESpace([Σ, V, Q]) + Y = MultiFieldFESpace([Σ, V, Q]) + + a((S, u, g), (T, v, e)) = ∫( S ⊙ T + u ⋅ (∇⋅T) + g * (χ ⊙ T) + + v ⋅ (∇⋅S) + e * (χ ⊙ S) )dΩ + l((T, v, e)) = ∫( -(fex ⋅ v) )dΩ + + solve(AffineFEOperator(a, l, X, Y)) +end + +S1, u1, γ1 = solve_weak_symmetry(1) +S2, u2, γ2 = solve_weak_symmetry(2) + +# ## Errors +# +# $O(h^k)$ in all three unknowns, so the $k=2$ element is one order better +# everywhere. + +for (k, Sh, uh, gh) in ((1, S1, u1, γ1), (2, S2, u2, γ2)) + eS = Sh - σex + eu = uh - uex + eg = gh - γex + println("k = ", k, + " |sigma| = ", sqrt(sum( ∫( eS ⊙ eS )dΩ )), + " |u| = ", sqrt(sum( ∫( eu ⋅ eu )dΩ )), + " |gamma| = ", sqrt(sum( ∫( eg * eg )dΩ ))) +end + +# ## How symmetric is it? +# +# Not pointwise, which is the price of the method: $\chi : \sigma_h$ is zero only +# against the multiplier space. Its $L^2$ norm converges to zero at the rate of +# the method, unlike for `aw_c`, where it is zero to round-off. + +for (k, Sh) in ((1, S1), (2, S2)) + asym = χ ⊙ Sh + println("k = ", k, " |skew part| = ", sqrt(sum( ∫( asym * asym )dΩ ))) +end + +# ## Visualisation + +mkpath("output_path") +writevtk(Ω, "output_path/weak_symmetry_elasticity", + cellfields=["sigma" => S2, "u" => u2, "gamma" => γ2]) From 1ccec7065ec47335cb31a2af3027f2276fd32a5d Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 16 Sep 2026 17:22:32 +1000 Subject: [PATCH 4/4] Activate new tests --- .github/workflows/ci.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bd38ff3..8897635 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -33,6 +33,10 @@ jobs: files: "stokes_blocks.jl poisson_amr.jl poisson_unfitted.jl" - name: Other files: "validation.jl validation_DrWatson.jl interpolation_fe.jl poisson_dev_fe.jl geometry_dev.jl arrays_dev.jl reffes_dev.jl tensor_values_dev.jl" + - name: Plates + files: "hermite_beam.jl morley_biharmonic.jl argyris_biharmonic.jl hhj_plate.jl laminate_plate.jl" + - name: Mixed + files: "arnold_winther_elasticity.jl weak_symmetry_elasticity.jl mtw_darcy_stokes.jl regge_metric.jl gls_stokes.jl" steps: - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v3