From af69c7ff51f4dcf248e3826bd8c306fae70110d5 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Sat, 29 Jan 2022 16:08:30 +0200 Subject: [PATCH] Add support for facet orientation in FacetBasis (#865) * Add workaround for default file format change in meshio 5.3.0 * Add skfem.generic_utils.OrientedBoundary * Add Mesh.facets_around for constructing OrientedBoundary * Rename BoundaryFacetBasis to FacetBasis and alias vice versa * Add example meshes with oriented interfaces * Support giving lambda AND index sets in with_boundaries * Fix the broken subdomain assembly with tags * Fix interpolation with multiple components * Allow specifying normal in Mesh.facets_satisfying * Support loading OrientedBoundary in Mesh.load from Gmsh format --- README.md | 34 +- docs/advanced.rst | 2 +- docs/api.rst | 8 +- docs/examples/meshes/cuubat.msh | 2473 +++++++++++++++++ docs/examples/meshes/interface.msh | 432 +++ docs/examples/meshes/internal.msh | 686 +++++ skfem/assembly/__init__.py | 2 +- skfem/assembly/basis/__init__.py | 6 +- skfem/assembly/basis/abstract_basis.py | 61 +- skfem/assembly/basis/cell_basis.py | 34 +- ...boundary_facet_basis.py => facet_basis.py} | 69 +- skfem/assembly/basis/interior_facet_basis.py | 38 +- skfem/assembly/basis/mortar_facet_basis.py | 32 +- skfem/element/__init__.py | 2 +- skfem/generic_utils.py | 17 + skfem/io/meshio.py | 75 +- skfem/mesh/mesh.py | 129 +- skfem/visuals/matplotlib.py | 9 +- tests/test_basis.py | 176 +- tests/test_mesh.py | 18 +- 20 files changed, 4119 insertions(+), 184 deletions(-) create mode 100644 docs/examples/meshes/cuubat.msh create mode 100644 docs/examples/meshes/interface.msh create mode 100644 docs/examples/meshes/internal.msh rename skfem/assembly/basis/{boundary_facet_basis.py => facet_basis.py} (83%) diff --git a/README.md b/README.md index ba58ec2a5..b9b8eeb4b 100644 --- a/README.md +++ b/README.md @@ -216,20 +216,44 @@ with respect to documented and/or tested features. ### Unreleased - Changed: `DiscreteField` is now a subclass of `ndarray` instead of - `NamedTuple` + `NamedTuple` and, consequently, the components of `DiscreteField` cannot no + more be indexed inside forms like `u[1]` (use `u.grad` instead) - Changed: Writing `w['u']` and `w.u` inside the form definition is now equivalent (previously `w.u == w['u'].value`) -- Changed: `Mesh.draw` now uses `matplotlib` by default, call - `Mesh.draw("vedo")` to keep using `vedo` +- Changed: `Mesh.draw` now uses `matplotlib` by default, calling + `Mesh.draw("vedo")` uses `vedo` - Changed: `skfem.visuals.matplotlib` now uses `jet` as the default colormap +- Changed: `BoundaryFacetBasis` is now an alias of `FacetBasis` instead of + other way around - Deprecated: `DiscreteField.value` remains for backwards-compatibility but is now deprecated and can be dropped - Added: `Mesh.plot`, a wrapper to `skfem.visuals.*.plot` - Added: `Basis.plot`, a wrapper to `skfem.visuals.*.plot` - Added: `Basis.refinterp` now supports vectorial fields -- Added: `skfem.visuals.matpltlib.plot` now has a basic quiver plot for vector +- Added: `skfem.visuals.matplotlib.plot` now has a basic quiver plot for vector fields -- Fixed: Improvements to backwards compatibility in `asm`/`assemble` kwargs +- Added: `Mesh.facets_around` which constructs a set of facets around a + subdomain +- Added: `Mesh.load` now tries loading the orientation of boundaries and + interfaces +- Added: `OrientedBoundary` which is a subclass of `ndarray` for facet index + arrays with the orientation information (0 or 1 per facet) available as + `OrientedBoundary.ori` +- Added: `FacetBasis` will use the facet orientations (if present) to calculate + traces and normal vectors +- Added: `skfem.visuals.matplotlib.draw` will visualize the orientations if + `boundaries=True` is given +- Added: `Mesh.facets_satisfying` allows specifying the keyword argument + `normal` for orienting the resulting interface +- Added: `FacetBasis` constructor now has the keyword argument `side` which + allows changing the side of the facet used to calculate the basis function + values and gradients +- Fixed: Improvements to backwards compatibility in `asm`/`assemble` keyword + arguments +- Fixed: Save format issue with meshio 5.3.0+ +- Fixed: `CellBasis` did not properly support `elements` argument +- Fixed: `Basis.interpolate` did not properly interpolate all components of + `ElementComposite` ### [5.2.0] - 2021-12-27 diff --git a/docs/advanced.rst b/docs/advanced.rst index 862c10a91..e9cf39b6b 100644 --- a/docs/advanced.rst +++ b/docs/advanced.rst @@ -65,7 +65,7 @@ This can be written as In addition, forms can depend on the local mesh parameter ``w.h`` or other finite element functions (see :ref:`predefined`). Moreover, boundary forms -assembled using :class:`~skfem.assembly.BoundaryFacetBasis` can depend on the +assembled using :class:`~skfem.assembly.FacetBasis` can depend on the outward normal vector ``w.n``. One example is the form .. math:: diff --git a/docs/api.rst b/docs/api.rst index b8c334771..3159b7078 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -90,12 +90,12 @@ Class: CellBasis :members: __init__, interpolate, project -Class: BoundaryFacetBasis -************************* - -.. autoclass:: skfem.assembly.FacetBasis +Class: FacetBasis +***************** .. autoclass:: skfem.assembly.BoundaryFacetBasis + +.. autoclass:: skfem.assembly.FacetBasis :members: __init__ Class: InteriorFacetBasis diff --git a/docs/examples/meshes/cuubat.msh b/docs/examples/meshes/cuubat.msh new file mode 100644 index 000000000..789bf0cdd --- /dev/null +++ b/docs/examples/meshes/cuubat.msh @@ -0,0 +1,2473 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +3 +2 26 "interface" +2 27 "boundary" +3 25 "both" +$EndPhysicalNames +$Entities +12 20 11 2 +1 0 0 1 0 +2 0 0 0 0 +3 0 1 1 0 +4 0 1 0 0 +5 1 0 0 0 +6 1 0 1 0 +7 1 1 1 0 +8 1 1 0 0 +9 2 0 1 0 +10 2 0 0 0 +11 2 1 1 0 +12 2 1 0 0 +1 -1e-07 -1e-07 -9.999999994736442e-08 1e-07 1e-07 1.0000001 0 2 2 -1 +2 -1e-07 -9.999999994736442e-08 0.9999999000000001 1e-07 1.0000001 1.0000001 0 2 1 -3 +3 -1e-07 0.9999999000000001 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 2 4 -3 +4 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 1.0000001 1e-07 0 2 2 -4 +5 0.9999999000000001 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 2 5 -6 +6 0.9999999000000001 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 6 -7 +7 0.9999999000000001 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 2 8 -7 +8 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 2 5 -8 +9 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 2 -5 +10 -9.999999994736442e-08 -1e-07 0.9999999000000001 1.0000001 1e-07 1.0000001 0 2 1 -6 +11 -9.999999994736442e-08 0.9999999000000001 -1e-07 1.0000001 1.0000001 1e-07 0 2 4 -8 +12 -9.999999994736442e-08 0.9999999000000001 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 3 -7 +13 1.9999999 -1e-07 -9.999999994736442e-08 2.0000001 1e-07 1.0000001 0 2 10 -9 +14 1.9999999 -9.999999994736442e-08 0.9999999000000001 2.0000001 1.0000001 1.0000001 0 2 9 -11 +15 1.9999999 0.9999999000000001 -9.999999994736442e-08 2.0000001 1.0000001 1.0000001 0 2 12 -11 +16 1.9999999 -9.999999994736442e-08 -1e-07 2.0000001 1.0000001 1e-07 0 2 10 -12 +17 0.9999999000000001 -1e-07 -1e-07 2.0000001 1e-07 1e-07 0 2 5 -10 +18 0.9999999000000001 -1e-07 0.9999999000000001 2.0000001 1e-07 1.0000001 0 2 6 -9 +19 0.9999999000000001 0.9999999000000001 -1e-07 2.0000001 1.0000001 1e-07 0 2 8 -12 +20 0.9999999000000001 0.9999999000000001 0.9999999000000001 2.0000001 1.0000001 1.0000001 0 2 7 -11 +1 -1e-07 -9.999999994736442e-08 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 1 27 4 1 2 -3 -4 +2 0.9999999000000001 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 1 26 4 5 6 -7 -8 +3 -9.999999994736442e-08 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 4 9 5 -10 -1 +4 -9.999999994736442e-08 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 11 7 -12 -3 +5 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 4 4 11 -8 -9 +6 -9.999999994736442e-08 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 4 2 12 -6 -10 +7 1.9999999 -9.999999994736442e-08 -9.999999994736442e-08 2.0000001 1.0000001 1.0000001 0 4 13 14 -15 -16 +8 0.9999999000000001 -1e-07 -9.999999994736442e-08 2.0000001 1e-07 1.0000001 0 4 17 13 -18 -5 +9 0.9999999000000001 0.9999999000000001 -9.999999994736442e-08 2.0000001 1.0000001 1.0000001 0 4 19 15 -20 -7 +10 0.9999999000000001 -9.999999994736442e-08 -1e-07 2.0000001 1.0000001 1e-07 0 4 8 19 -16 -17 +11 0.9999999000000001 -9.999999994736442e-08 0.9999999000000001 2.0000001 1.0000001 1.0000001 0 4 6 20 -14 -18 +1 -9.999999994736442e-08 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 1 25 6 1 2 3 4 5 6 +2 0.9999999000000001 -9.999999994736442e-08 -9.999999994736442e-08 2.0000001 1.0000001 1.0000001 1 25 6 2 7 8 9 10 11 +$EndEntities +$Nodes +45 419 1 419 +0 1 0 1 +1 +0 0 1 +0 2 0 1 +2 +0 0 0 +0 3 0 1 +3 +0 1 1 +0 4 0 1 +4 +0 1 0 +0 5 0 1 +5 +1 0 0 +0 6 0 1 +6 +1 0 1 +0 7 0 1 +7 +1 1 1 +0 8 0 1 +8 +1 1 0 +0 9 0 1 +9 +2 0 1 +0 10 0 1 +10 +2 0 0 +0 11 0 1 +11 +2 1 1 +0 12 0 1 +12 +2 1 0 +1 1 0 4 +13 +14 +15 +16 +0 0 0.2000000000000003 +0 0 0.4000000000000009 +0 0 0.6000000000000009 +0 0 0.8000000000000003 +1 2 0 4 +17 +18 +19 +20 +0 0.2000000000000003 1 +0 0.4000000000000009 1 +0 0.6000000000000009 1 +0 0.8000000000000003 1 +1 3 0 4 +21 +22 +23 +24 +0 1 0.2000000000000003 +0 1 0.4000000000000009 +0 1 0.6000000000000009 +0 1 0.8000000000000003 +1 4 0 4 +25 +26 +27 +28 +0 0.2000000000000003 0 +0 0.4000000000000009 0 +0 0.6000000000000009 0 +0 0.8000000000000003 0 +1 5 0 4 +29 +30 +31 +32 +1 0 0.2000000000000003 +1 0 0.4000000000000009 +1 0 0.6000000000000009 +1 0 0.8000000000000003 +1 6 0 4 +33 +34 +35 +36 +1 0.2000000000000003 1 +1 0.4000000000000009 1 +1 0.6000000000000009 1 +1 0.8000000000000003 1 +1 7 0 4 +37 +38 +39 +40 +1 1 0.2000000000000003 +1 1 0.4000000000000009 +1 1 0.6000000000000009 +1 1 0.8000000000000003 +1 8 0 4 +41 +42 +43 +44 +1 0.2000000000000003 0 +1 0.4000000000000009 0 +1 0.6000000000000009 0 +1 0.8000000000000003 0 +1 9 0 4 +45 +46 +47 +48 +0.2000000000000003 0 0 +0.4000000000000009 0 0 +0.6000000000000009 0 0 +0.8000000000000003 0 0 +1 10 0 4 +49 +50 +51 +52 +0.2000000000000003 0 1 +0.4000000000000009 0 1 +0.6000000000000009 0 1 +0.8000000000000003 0 1 +1 11 0 4 +53 +54 +55 +56 +0.2000000000000003 1 0 +0.4000000000000009 1 0 +0.6000000000000009 1 0 +0.8000000000000003 1 0 +1 12 0 4 +57 +58 +59 +60 +0.2000000000000003 1 1 +0.4000000000000009 1 1 +0.6000000000000009 1 1 +0.8000000000000003 1 1 +1 13 0 4 +61 +62 +63 +64 +2 0 0.2000000000000003 +2 0 0.4000000000000009 +2 0 0.6000000000000009 +2 0 0.8000000000000003 +1 14 0 4 +65 +66 +67 +68 +2 0.2000000000000003 1 +2 0.4000000000000009 1 +2 0.6000000000000009 1 +2 0.8000000000000003 1 +1 15 0 4 +69 +70 +71 +72 +2 1 0.2000000000000003 +2 1 0.4000000000000009 +2 1 0.6000000000000009 +2 1 0.8000000000000003 +1 16 0 4 +73 +74 +75 +76 +2 0.2000000000000003 0 +2 0.4000000000000009 0 +2 0.6000000000000009 0 +2 0.8000000000000003 0 +1 17 0 4 +77 +78 +79 +80 +1.2 0 0 +1.400000000000001 0 0 +1.600000000000001 0 0 +1.8 0 0 +1 18 0 4 +81 +82 +83 +84 +1.2 0 1 +1.400000000000001 0 1 +1.600000000000001 0 1 +1.8 0 1 +1 19 0 4 +85 +86 +87 +88 +1.2 1 0 +1.400000000000001 1 0 +1.600000000000001 1 0 +1.8 1 0 +1 20 0 4 +89 +90 +91 +92 +1.2 1 1 +1.400000000000001 1 1 +1.600000000000001 1 1 +1.8 1 1 +2 1 0 24 +93 +94 +95 +96 +97 +98 +99 +100 +101 +102 +103 +104 +105 +106 +107 +108 +109 +110 +111 +112 +113 +114 +115 +116 +0 0.8267949192431125 0.5000000000000009 +0 0.1335085296108589 0.4919480156549869 +0 0.4903320945051429 0.1627551561799211 +0 0.4886650857412674 0.8362933879517429 +0 0.8306455909888337 0.3008474738561762 +0 0.6705807092817124 0.3964399829492918 +0 0.6727576366574796 0.6030990510214808 +0 0.4803847577293379 0.500000000000001 +0 0.3001135766047719 0.8306178565232263 +0 0.1880910040222625 0.3152349128312829 +0 0.8303769882505059 0.6995291129326942 +0 0.3057723032641354 0.1616307679810649 +0 0.6838607148057155 0.1881003643263491 +0 0.1745203995428645 0.6841099050678019 +0 0.6856863617897432 0.817662137456531 +0 0.3853794790740478 0.6683776366978199 +0 0.3888753081903478 0.3199030520313019 +0 0.3176836966063226 0.4765884864877661 +0 0.8535898384862238 0.1464101615137761 +0 0.1464101615137758 0.8535898384862242 +0 0.1464101615137758 0.1464101615137758 +0 0.8535898384862249 0.8535898384862249 +0 0.5417981355829327 0.6896969712338845 +0 0.5418767703780959 0.3102497757852787 +2 2 0 24 +117 +118 +119 +120 +121 +122 +123 +124 +125 +126 +127 +128 +129 +130 +131 +132 +133 +134 +135 +136 +137 +138 +139 +140 +1 0.5000000000000009 0.8267949192431124 +1 0.491948015654987 0.1335085296108589 +1 0.1627551561799211 0.4903320945051428 +1 0.8362933879517429 0.4886650857412675 +1 0.3008474738561758 0.8306455909888335 +1 0.3964399829492918 0.6705807092817122 +1 0.6030990510214806 0.6727576366574792 +1 0.5000000000000011 0.4803847577293376 +1 0.315234912831283 0.1880910040222625 +1 0.8306178565232264 0.3001135766047722 +1 0.6995291129326939 0.8303769882505059 +1 0.161630767981065 0.3057723032641353 +1 0.6841099050678021 0.1745203995428646 +1 0.1881003643263491 0.6838607148057156 +1 0.8176621374565309 0.6856863617897432 +1 0.6683776366978198 0.3853794790740478 +1 0.3199030520313019 0.3888753081903477 +1 0.476588486487766 0.3176836966063225 +1 0.8535898384862239 0.146410161513776 +1 0.146410161513776 0.146410161513776 +1 0.1464101615137758 0.8535898384862242 +1 0.8535898384862249 0.8535898384862249 +1 0.6896969712338844 0.5417981355829328 +1 0.3102497757852788 0.5418767703780956 +2 3 0 24 +141 +142 +143 +144 +145 +146 +147 +148 +149 +150 +151 +152 +153 +154 +155 +156 +157 +158 +159 +160 +161 +162 +163 +164 +0.8267949192431123 0 0.5000000000000007 +0.1335085296108589 0 0.491948015654987 +0.4903320945051427 0 0.162755156179921 +0.4886650857412677 0 0.8362933879517431 +0.8306455909888335 0 0.3008474738561758 +0.670580709281712 0 0.3964399829492917 +0.6727576366574792 0 0.6030990510214806 +0.4803847577293376 0 0.5000000000000011 +0.1880910040222625 0 0.3152349128312829 +0.8303769882505057 0 0.6995291129326938 +0.3001135766047722 0 0.8306178565232264 +0.3057723032641353 0 0.161630767981065 +0.1745203995428646 0 0.684109905067802 +0.6838607148057154 0 0.188100364326349 +0.6856863617897432 0 0.8176621374565309 +0.3853794790740479 0 0.6683776366978199 +0.3888753081903476 0 0.3199030520313019 +0.3176836966063225 0 0.4765884864877659 +0.146410161513776 0 0.8535898384862239 +0.8535898384862242 0 0.1464101615137758 +0.146410161513776 0 0.146410161513776 +0.8535898384862247 0 0.8535898384862247 +0.5417981355829328 0 0.6896969712338845 +0.5418767703780953 0 0.3102497757852787 +2 4 0 24 +165 +166 +167 +168 +169 +170 +171 +172 +173 +174 +175 +176 +177 +178 +179 +180 +181 +182 +183 +184 +185 +186 +187 +188 +0.5000000000000009 1 0.8267949192431125 +0.4939184515869166 1 0.1323760430703402 +0.1637066120482572 1 0.4886650857412678 +0.8362933879517429 1 0.4886650857412674 +0.3008474738561762 1 0.8306455909888337 +0.3964488584460981 1 0.6705676034825185 +0.6031005302709485 1 0.6727554523576138 +0.500000000000001 1 0.4803847577293379 +0.1693821434767736 1 0.3001135766047722 +0.8306178565232263 1 0.3001135766047719 +0.6995291771094898 1 0.8303768934853097 +0.3157845613177561 1 0.1865278143474288 +0.1882451620900605 1 0.6836094648973968 +0.6859226860940238 1 0.1724846426120655 +0.8176622949490069 1 0.6856861292320839 +0.3214614339239139 1 0.3873807148748087 +0.6687225910033419 1 0.3893144187516059 +0.4768920018051231 1 0.3175842984916424 +0.1464101615137761 1 0.8535898384862238 +0.8535898384862242 1 0.1464101615137758 +0.1464101615137759 1 0.1464101615137759 +0.8535898384862249 1 0.8535898384862249 +0.310303028766116 1 0.5417981355829332 +0.6896969712338845 1 0.5417981355829327 +2 5 0 24 +189 +190 +191 +192 +193 +194 +195 +196 +197 +198 +199 +200 +201 +202 +203 +204 +205 +206 +207 +208 +209 +210 +211 +212 +0.8267949192431125 0.5000000000000009 0 +0.1335085296108589 0.4919480156549869 0 +0.4903320945051429 0.1627551561799211 0 +0.4886650857412674 0.8362933879517429 0 +0.8306455909888337 0.3008474738561762 0 +0.6705807092817124 0.3964399829492918 0 +0.6727576366574796 0.6030990510214808 0 +0.4803847577293379 0.500000000000001 0 +0.1880910040222625 0.3152349128312829 0 +0.3001135766047719 0.8306178565232263 0 +0.8303769882505059 0.6995291129326942 0 +0.3057723032641354 0.1616307679810649 0 +0.6838607148057155 0.1881003643263491 0 +0.1745203995428645 0.6841099050678019 0 +0.6856863617897432 0.817662137456531 0 +0.3853794790740478 0.6683776366978199 0 +0.3888753081903478 0.3199030520313019 0 +0.3176836966063226 0.4765884864877661 0 +0.8535898384862238 0.1464101615137761 0 +0.1464101615137758 0.8535898384862242 0 +0.1464101615137758 0.1464101615137758 0 +0.8535898384862249 0.8535898384862249 0 +0.5417981355829327 0.6896969712338845 0 +0.5418767703780959 0.3102497757852787 0 +2 6 0 25 +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 +0.4983653964774457 0.172261341973122 1 +0.5118616330494621 0.8419191893773027 1 +0.1627551561799213 0.5096679054948585 1 +0.8363011050363136 0.5080333019723035 1 +0.3007481903574787 0.1692970876564672 1 +0.3961510017790835 0.3292524473618773 1 +0.6002089867631721 0.3244658182136438 1 +0.5000000000000012 0.5196152422706624 1 +0.3180797153162082 0.8083986701189039 1 +0.8148028884944708 0.3106709346301967 1 +0.69992182039465 0.8517369372914835 1 +0.1880448979450769 0.3161072616641265 1 +0.6944945183702189 0.1589680427217897 1 +0.161969434943556 0.6938098008002902 1 +0.8086055710882908 0.7054234997311373 1 +0.3204336302725384 0.6104699881772513 1 +0.6641620538129379 0.5994423454100184 1 +0.469547699410264 0.6732050807568873 1 +0.8535898384862239 0.146410161513776 1 +0.146410161513776 0.146410161513776 1 +0.1464101615137758 0.8535898384862242 1 +0.8535898384862248 0.8535898384862248 1 +0.3102497757852793 0.4581232296219057 1 +0.6880178165217438 0.4566452774346555 1 +0.6092224884805854 0.7272257338607353 1 +2 7 0 25 +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 +2 0.4983653964774457 0.1722613419731219 +2 0.5118616330494621 0.8419191893773028 +2 0.1627551561799212 0.5096679054948585 +2 0.8363011050363134 0.5080333019723035 +2 0.3007481903574787 0.1692970876564672 +2 0.3961510017790833 0.3292524473618773 +2 0.6002089867631721 0.3244658182136438 +2 0.5000000000000011 0.5196152422706624 +2 0.318079715316208 0.8083986701189039 +2 0.8148028884944708 0.3106709346301967 +2 0.6999218203946499 0.8517369372914836 +2 0.1880448979450769 0.3161072616641265 +2 0.6944945183702189 0.1589680427217897 +2 0.161969434943556 0.6938098008002902 +2 0.8086055710882907 0.7054234997311373 +2 0.3204336302725382 0.6104699881772513 +2 0.6641620538129378 0.5994423454100184 +2 0.4695476994102636 0.6732050807568875 +2 0.8535898384862239 0.146410161513776 +2 0.146410161513776 0.146410161513776 +2 0.1464101615137758 0.8535898384862242 +2 0.8535898384862248 0.8535898384862248 +2 0.3102497757852793 0.4581232296219056 +2 0.6880178165217437 0.4566452774346555 +2 0.609222488480585 0.7272257338607353 +2 8 0 24 +263 +264 +265 +266 +267 +268 +269 +270 +271 +272 +273 +274 +275 +276 +277 +278 +279 +280 +281 +282 +283 +284 +285 +286 +1.173205080756887 0 0.5000000000000009 +1.866491470389141 0 0.4919480156549869 +1.509667905494859 0 0.1627551561799213 +1.511334914258734 0 0.8362933879517429 +1.169354409011166 0 0.3008474738561762 +1.329419290718288 0 0.3964399829492921 +1.327242363342521 0 0.6030990510214808 +1.519615242270662 0 0.500000000000001 +1.811908995977738 0 0.315234912831283 +1.169623011749494 0 0.6995291129326938 +1.699886423395229 0 0.8306178565232268 +1.316139285194285 0 0.1881003643263493 +1.314313638210258 0 0.8176621374565307 +1.825479600457136 0 0.6841099050678017 +1.694227696735866 0 0.1616307679810649 +1.614620520925953 0 0.6683776366978199 +1.611124691809653 0 0.3199030520313021 +1.682316303393677 0 0.4765884864877665 +1.853589838486224 0 0.1464101615137758 +1.146410161513776 0 0.146410161513776 +1.853589838486225 0 0.8535898384862247 +1.146410161513775 0 0.8535898384862247 +1.458201864417069 0 0.689696971233884 +1.458123229621906 0 0.3102497757852794 +2 9 0 24 +287 +288 +289 +290 +291 +292 +293 +294 +295 +296 +297 +298 +299 +300 +301 +302 +303 +304 +305 +306 +307 +308 +309 +310 +1.827738658026878 1 0.4983653964774455 +1.133508529610859 1 0.491948015654987 +1.490332094505143 1 0.1627551561799211 +1.490299689263823 1 0.8353496491679775 +1.830702912343533 1 0.3007481903574784 +1.670747552638123 1 0.396151001779083 +1.675521075987163 1 0.6002001112663657 +1.480384757729337 1 0.5000000000000009 +1.188091004022263 1 0.3152349128312829 +1.300184493389556 1 0.8305769126984412 +1.689149601149577 1 0.8146994615203915 +1.305772303264135 1 0.161630767981065 +1.174524091587054 1 0.6841077734650955 +1.840999862274974 1 0.6944758012917381 +1.683892738335874 1 0.1880448979450768 +1.385664347799302 1 0.6682131676626103 +1.388875308190348 1 0.3199030520313018 +1.317696215728973 1 0.4765812585689336 +1.146410161513776 1 0.8535898384862239 +1.853589838486224 1 0.1464101615137758 +1.146410161513776 1 0.146410161513776 +1.853589838486225 1 0.8535898384862247 +1.541798135582933 1 0.6896969712338844 +1.541876770378095 1 0.3102497757852788 +2 10 0 24 +311 +312 +313 +314 +315 +316 +317 +318 +319 +320 +321 +322 +323 +324 +325 +326 +327 +328 +329 +330 +331 +332 +333 +334 +1.827738658026878 0.4983653964774456 0 +1.133508529610859 0.4919480156549869 0 +1.490332094505143 0.1627551561799211 0 +1.490299689263823 0.8353496491679776 0 +1.830702912343533 0.3007481903574787 0 +1.670747552638123 0.3961510017790832 0 +1.675521075987163 0.6002001112663659 0 +1.480384757729338 0.5000000000000011 0 +1.188091004022263 0.3152349128312829 0 +1.689149601149577 0.8146994615203915 0 +1.300184493389556 0.8305769126984412 0 +1.305772303264135 0.1616307679810649 0 +1.174524091587054 0.6841077734650954 0 +1.683892738335874 0.1880448979450768 0 +1.840999862274974 0.6944758012917381 0 +1.385664347799302 0.6682131676626104 0 +1.388875308190348 0.3199030520313019 0 +1.317696215728973 0.4765812585689336 0 +1.853589838486224 0.1464101615137761 0 +1.146410161513776 0.8535898384862244 0 +1.146410161513776 0.1464101615137758 0 +1.853589838486225 0.853589838486225 0 +1.541798135582933 0.6896969712338846 0 +1.541876770378096 0.3102497757852787 0 +2 11 0 24 +335 +336 +337 +338 +339 +340 +341 +342 +343 +344 +345 +346 +347 +348 +349 +350 +351 +352 +353 +354 +355 +356 +357 +358 +1.827738658026878 0.4983653964774455 1 +1.133508529610859 0.4919480156549869 1 +1.490332094505143 0.1627551561799211 1 +1.490299689263823 0.8353496491679776 1 +1.830702912343533 0.3007481903574788 1 +1.670747552638123 0.3961510017790831 1 +1.675521075987163 0.6002001112663659 1 +1.480384757729337 0.5000000000000011 1 +1.188091004022263 0.3152349128312829 1 +1.689149601149577 0.8146994615203915 1 +1.300184493389556 0.8305769126984412 1 +1.305772303264135 0.1616307679810649 1 +1.683892738335874 0.1880448979450768 1 +1.174524091587054 0.6841077734650954 1 +1.840999862274974 0.6944758012917381 1 +1.385664347799302 0.6682131676626103 1 +1.388875308190348 0.3199030520313019 1 +1.317696215728973 0.4765812585689335 1 +1.853589838486224 0.1464101615137761 1 +1.146410161513776 0.8535898384862244 1 +1.853589838486225 0.8535898384862249 1 +1.146410161513776 0.1464101615137758 1 +1.541798135582933 0.6896969712338845 1 +1.541876770378096 0.3102497757852788 1 +3 1 0 30 +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 +0.5 0.5 0.4999999999999999 +0.3137161618329791 0.3137161618329792 0.3137161618329791 +0.68716108135786 0.6850327481252375 0.6885882840511444 +0.3119844512131699 0.6855505714011596 0.3140971213318855 +0.6419384953854816 0.3000211161636266 0.3000211161636269 +0.312659123099453 0.6876671548417834 0.68570428033854 +0.3143725187296867 0.3143725187296866 0.6879940089604385 +0.6856156542067544 0.6864906145802389 0.3121124389250448 +0.6870349332794939 0.3128449153779588 0.6876646984044233 +0.2496475171256845 0.499177301646953 0.4717494239768313 +0.7495718051243533 0.4475098633020375 0.4792775283178745 +0.4693290054628529 0.5004267570388961 0.2508184802477197 +0.5015822188356445 0.7515446535274326 0.5295661858796084 +0.4690507314701258 0.2508293970712898 0.500694114208053 +0.4705396416389634 0.5021297695005823 0.7503012058398677 +0.2167079917035499 0.2167079917035499 0.5011964755387098 +0.224724075223921 0.4993283333688674 0.2110419722589371 +0.4980756649900309 0.7890714622930484 0.2255433708839236 +0.2206044085812091 0.5011615202940687 0.7794730912381842 +0.778520289889941 0.778520289889941 0.5003313602421918 +0.4998140187670144 0.7762460508250575 0.7817392021850544 +0.7802737094476384 0.2212218141795319 0.4644107389728864 +0.222323847984227 0.7809701555553022 0.4995266378866139 +0.5003964799876288 0.2193011318360893 0.777609472145723 +0.7777874731988524 0.4986776632869084 0.7809856898418527 +0.7782391794437367 0.4826335509748605 0.2321146530088503 +0.4696604487881277 0.2090028699660801 0.2090028699660802 +0.8085019154388607 0.1952134768551254 0.2221569063338013 +0.2231590239095343 0.1847211483465517 0.1847211483465517 +0.7772287115298642 0.8134455096342376 0.8167431129612636 +3 2 0 31 +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 +1.5 0.4999999999999999 0.5 +1.314209177684468 0.3120432064403601 0.314209177684468 +1.688410772828302 0.6875900252735611 0.6852171695321694 +1.313946010553994 0.6863958723835584 0.3138853775822013 +1.686327107505047 0.3144584552580744 0.6877729388705214 +1.641779330384691 0.6999810943264821 0.299991814110578 +1.687389679131306 0.3143602514708858 0.3126103208686937 +1.312295900604692 0.312646888534806 0.6855830402693648 +1.312069198486443 0.6858642882546133 0.6858642882546134 +1.752311053429209 0.5222783513688078 0.4508979885822303 +1.472701096031907 0.2485442930405478 0.5008212808238963 +1.249642261560278 0.4705252601187587 0.4993112680118318 +1.499550786661745 0.4702977897775137 0.2496328651275853 +1.529026918739945 0.5023538771055848 0.7515234978162464 +1.499680994128748 0.7503722429744359 0.5301796381198702 +1.216764235931495 0.4990669829023401 0.2167642359314946 +1.77480240512919 0.7705453725216593 0.4817934796621313 +1.774377995473458 0.2117344372048864 0.5004663878189919 +1.224898548879484 0.7889885268237906 0.499577272443817 +1.500976344743215 0.2205522314216368 0.2204362386760998 +1.499218310021813 0.220820819500757 0.7793066686055991 +1.21879888705506 0.2223225376041431 0.4995558291590242 +1.781920402528074 0.5011150101414229 0.7768329051629101 +1.77656814085874 0.5059106421428522 0.2234318591412597 +1.499979262144692 0.7784082057864308 0.7784082057864308 +1.226787357329053 0.4984679970811196 0.7901579956508561 +1.469661685351805 0.7912327873149936 0.2092960759702209 +1.807539232584233 0.7771785026404542 0.1924577898127455 +1.815458977712783 0.7765518489238763 0.8161353523670102 +1.222918868925507 0.8164033191582545 0.1872356087715577 +1.183426734106361 0.1867317481649642 0.2229669924253743 +$EndNodes +$Elements +4 1523 1 1523 +2 1 2 66 +1 16 1 112 +2 1 17 112 +3 2 13 113 +4 25 2 113 +5 20 3 114 +6 3 24 114 +7 21 4 111 +8 4 28 111 +9 13 14 102 +10 13 102 113 +11 14 15 94 +12 14 94 102 +13 15 16 106 +14 94 15 106 +15 106 16 112 +16 17 18 101 +17 17 101 112 +18 18 19 96 +19 18 96 101 +20 19 20 107 +21 96 19 107 +22 107 20 114 +23 22 21 97 +24 97 21 111 +25 23 22 93 +26 93 22 97 +27 24 23 103 +28 23 93 103 +29 24 103 114 +30 26 25 104 +31 104 25 113 +32 27 26 95 +33 95 26 104 +34 28 27 105 +35 27 95 105 +36 28 105 111 +37 93 97 98 +38 93 98 99 +39 93 99 103 +40 102 94 110 +41 94 106 110 +42 95 104 109 +43 105 95 116 +44 95 109 116 +45 101 96 108 +46 96 107 115 +47 108 96 115 +48 98 97 105 +49 105 97 111 +50 99 98 100 +51 100 98 116 +52 98 105 116 +53 99 100 115 +54 103 99 107 +55 107 99 115 +56 108 100 110 +57 100 108 115 +58 100 109 110 +59 109 100 116 +60 106 101 108 +61 101 106 112 +62 104 102 109 +63 102 104 113 +64 109 102 110 +65 103 107 114 +66 106 108 110 +2 2 2 66 +67 5 136 29 +68 41 136 5 +69 32 137 6 +70 6 137 33 +71 36 138 7 +72 7 138 40 +73 37 135 8 +74 8 135 44 +75 29 128 30 +76 29 136 128 +77 30 119 31 +78 30 128 119 +79 31 130 32 +80 119 130 31 +81 130 137 32 +82 33 121 34 +83 33 137 121 +84 34 117 35 +85 34 121 117 +86 35 127 36 +87 117 127 35 +88 127 138 36 +89 38 126 37 +90 126 135 37 +91 39 120 38 +92 120 126 38 +93 40 131 39 +94 39 131 120 +95 40 138 131 +96 42 125 41 +97 125 136 41 +98 43 118 42 +99 118 125 42 +100 44 129 43 +101 43 129 118 +102 44 135 129 +103 121 122 117 +104 122 123 117 +105 123 127 117 +106 118 134 125 +107 129 134 118 +108 128 133 119 +109 119 140 130 +110 133 140 119 +111 120 132 126 +112 131 139 120 +113 120 139 132 +114 121 130 122 +115 121 137 130 +116 122 124 123 +117 122 140 124 +118 130 140 122 +119 124 139 123 +120 123 131 127 +121 123 139 131 +122 124 134 132 +123 132 139 124 +124 133 134 124 +125 124 140 133 +126 125 133 128 +127 128 136 125 +128 125 134 133 +129 126 132 129 +130 129 135 126 +131 131 138 127 +132 132 134 129 +3 1 4 690 +133 173 97 111 362 +134 217 151 232 365 +135 184 126 174 366 +136 233 226 114 364 +137 198 185 208 362 +138 111 208 185 362 +139 185 173 111 362 +140 159 112 232 365 +141 232 112 101 365 +142 159 232 151 365 +143 114 183 233 364 +144 184 135 126 366 +145 210 135 184 366 +146 169 233 183 364 +147 114 103 183 364 +148 199 135 210 366 +149 137 162 231 367 +150 225 231 162 367 +151 107 114 226 364 +152 231 121 137 367 +153 137 150 162 367 +154 176 185 198 362 +155 225 162 155 367 +156 101 224 232 365 +157 222 121 231 367 +158 177 183 103 364 +159 233 169 221 364 +160 137 130 150 367 +161 176 173 185 362 +162 114 107 103 364 +163 210 203 199 366 +164 224 217 232 365 +165 105 111 97 362 +166 225 222 231 367 +167 233 221 226 364 +168 162 150 155 367 +169 169 183 177 364 +170 137 121 130 367 +171 112 159 153 365 +172 129 135 199 366 +173 359 361 369 378 +174 359 369 366 378 +175 361 359 371 378 +176 361 359 369 383 +177 366 371 359 378 +178 359 361 373 383 +179 367 369 359 383 +180 359 373 367 383 +181 202 198 208 362 +182 101 112 106 365 +183 151 153 159 365 +184 126 135 129 366 +185 97 362 173 381 +186 174 366 126 378 +187 151 365 217 382 +188 184 174 178 366 +189 180 371 362 376 +190 362 371 370 376 +191 364 115 368 377 +192 370 371 366 376 +193 367 373 236 383 +194 150 367 130 380 +195 124 139 132 378 +196 367 372 163 382 +197 115 99 364 368 +198 226 107 364 377 +199 366 369 132 378 +200 219 236 367 373 +201 211 370 366 376 +202 108 368 365 374 +203 176 362 198 376 +204 163 367 147 372 +205 362 371 180 381 +206 112 153 106 365 +207 211 195 366 370 +208 366 370 359 371 +209 104 109 360 375 +210 224 365 101 377 +211 225 367 155 382 +212 361 373 229 379 +213 362 359 370 371 +214 177 364 103 381 +215 121 367 222 383 +216 124 132 369 378 +217 366 371 181 376 +218 169 221 364 379 +219 147 372 367 380 +220 104 360 109 387 +221 104 375 360 387 +222 365 368 108 377 +223 132 369 366 384 +224 157 372 360 374 +225 226 364 228 377 +226 133 380 384 386 +227 150 147 367 380 +228 369 370 366 384 +229 366 369 359 370 +230 229 373 361 383 +231 368 371 362 381 +232 229 373 230 379 +233 364 371 368 381 +234 205 360 370 375 +235 173 362 180 381 +236 228 226 221 364 +237 115 108 368 377 +238 173 176 180 362 +239 188 361 371 378 +240 188 179 361 378 +241 367 373 372 382 +242 187 371 364 381 +243 198 362 204 376 +244 146 372 147 380 +245 109 102 360 374 +246 149 157 360 374 +247 99 364 368 381 +248 107 115 364 377 +249 195 366 370 384 +250 133 380 369 384 +251 146 363 164 372 +252 131 139 361 378 +253 194 212 363 370 +254 182 371 180 376 +255 237 229 230 379 +256 101 365 108 377 +257 368 373 364 377 +258 365 373 368 377 +259 363 370 369 384 +260 181 371 366 378 +261 219 373 367 382 +262 229 236 373 383 +263 171 188 361 371 +264 205 197 360 375 +265 132 126 366 378 +266 116 362 368 375 +267 179 171 188 361 +268 156 365 151 382 +269 203 211 366 376 +270 187 170 364 371 +271 130 367 140 380 +272 155 367 163 382 +273 225 219 222 367 +274 372 373 365 382 +275 237 361 229 379 +276 359 372 367 373 +277 235 373 365 377 +278 146 363 372 380 +279 105 362 116 375 +280 222 367 236 383 +281 103 107 99 364 +282 139 123 131 361 +283 99 107 115 364 +284 195 370 194 384 +285 187 180 371 381 +286 170 169 364 379 +287 98 362 97 381 +288 164 372 363 385 +289 212 363 370 385 +290 98 116 362 368 +291 181 366 174 378 +292 359 362 368 371 +293 219 367 225 382 +294 180 362 176 376 +295 360 372 157 385 +296 205 370 360 385 +297 187 364 177 381 +298 218 217 365 382 +299 157 149 360 387 +300 205 360 197 387 +301 199 203 195 366 +302 235 365 224 377 +303 150 130 141 380 +304 359 368 364 371 +305 363 359 369 370 +306 147 155 150 367 +307 156 365 372 374 +308 99 103 364 381 +309 363 384 380 386 +310 195 199 366 384 +311 155 147 163 367 +312 129 366 199 384 +313 218 365 235 373 +314 203 211 195 366 +315 140 367 122 369 +316 97 98 105 362 +317 363 194 370 384 +318 217 224 218 365 +319 177 170 169 364 +320 109 368 360 375 +321 219 236 222 367 +322 228 364 221 379 +323 140 367 369 380 +324 204 370 362 375 +325 122 367 121 383 +326 109 360 102 387 +327 363 370 360 372 +328 105 98 116 362 +329 130 121 122 367 +330 187 170 177 364 +331 100 115 108 368 +332 224 235 218 365 +333 171 361 179 388 +334 141 130 119 380 +335 140 130 122 367 +336 214 175 223 379 +337 164 157 372 385 +338 359 364 368 373 +339 212 370 205 385 +340 367 372 369 380 +341 363 385 201 386 +342 359 368 365 373 +343 154 201 385 386 +344 360 375 197 387 +345 174 126 120 378 +346 362 370 204 376 +347 179 361 378 388 +348 229 220 230 373 +349 148 147 146 372 +350 104 197 375 387 +351 196 194 195 370 +352 174 120 168 378 +353 108 110 368 374 +354 154 385 363 386 +355 151 217 144 382 +356 149 374 360 387 +357 105 202 362 375 +358 198 202 204 362 +359 97 173 167 381 +360 211 204 370 376 +361 227 127 383 388 +362 122 369 367 383 +363 364 373 228 377 +364 101 106 108 365 +365 151 156 153 365 +366 157 158 372 374 +367 131 361 123 388 +368 191 143 154 385 +369 359 370 363 372 +370 178 203 366 376 +371 369 372 363 380 +372 220 229 236 373 +373 132 126 129 366 +374 359 360 370 372 +375 369 380 363 384 +376 182 181 371 376 +377 193 384 194 386 +378 360 374 102 387 +379 156 372 365 382 +380 127 131 123 388 +381 229 383 361 388 +382 359 369 367 372 +383 181 188 371 378 +384 98 368 362 381 +385 172 187 180 371 +386 228 221 230 379 +387 99 368 98 381 +388 144 217 213 382 +389 175 223 379 388 +390 359 365 372 373 +391 127 117 216 383 +392 368 372 365 374 +393 104 190 95 375 +394 214 165 175 379 +395 167 93 97 381 +396 178 174 181 366 +397 104 109 102 387 +398 365 373 218 382 +399 219 218 373 382 +400 229 227 383 388 +401 131 378 361 388 +402 216 227 127 383 +403 106 108 365 374 +404 156 153 365 374 +405 360 372 368 374 +406 132 134 369 384 +407 359 361 371 373 +408 371 373 361 379 +409 168 120 179 378 +410 204 362 202 375 +411 194 384 363 386 +412 156 163 372 382 +413 179 378 131 388 +414 219 220 236 373 +415 365 368 359 372 +416 115 100 99 368 +417 117 121 222 383 +418 182 181 172 371 +419 150 141 147 380 +420 205 370 206 375 +421 146 164 148 372 +422 364 371 170 379 +423 171 371 361 379 +424 129 132 366 384 +425 194 196 212 370 +426 359 363 369 372 +427 163 147 148 372 +428 211 196 195 370 +429 133 384 125 386 +430 222 216 117 383 +431 116 368 109 375 +432 213 225 155 382 +433 237 229 361 388 +434 226 215 107 377 +435 149 102 374 387 +436 93 177 103 381 +437 164 363 154 385 +438 201 363 212 385 +439 224 101 96 377 +440 179 168 39 120 +441 155 144 213 382 +442 178 366 181 376 +443 132 134 124 369 +444 364 373 371 379 +445 182 172 180 371 +446 169 165 221 379 +447 106 365 153 374 +448 201 194 363 386 +449 154 201 191 385 +450 190 105 95 375 +451 363 372 360 385 +452 360 370 363 385 +453 119 128 145 380 +454 167 177 93 381 +455 360 385 157 387 +456 205 385 360 387 +457 164 146 154 363 +458 201 194 212 363 +459 360 359 368 372 +460 221 165 214 379 +461 193 118 189 384 +462 94 142 149 374 +463 197 190 104 375 +464 203 166 192 376 +465 47 191 143 154 +466 148 164 157 372 +467 179 120 131 378 +468 212 196 205 370 +469 100 108 110 368 +470 359 371 364 373 +471 171 379 361 388 +472 117 34 222 121 +473 117 34 216 222 +474 145 141 119 380 +475 128 30 119 145 +476 155 51 225 213 +477 215 96 107 377 +478 51 155 144 213 +479 201 193 194 386 +480 235 228 373 377 +481 109 110 102 374 +482 149 158 157 374 +483 215 224 96 377 +484 181 188 172 371 +485 99 100 98 368 +486 223 59 175 214 +487 359 360 368 370 +488 145 30 119 141 +489 217 50 151 144 +490 220 219 218 373 +491 189 129 199 384 +492 39 131 179 120 +493 154 363 146 386 +494 193 125 118 384 +495 141 31 119 130 +496 361 379 237 388 +497 156 372 158 374 +498 97 173 22 167 +499 217 50 144 213 +500 191 200 152 385 +501 125 384 193 386 +502 31 141 150 130 +503 93 97 22 167 +504 103 93 23 177 +505 167 23 93 177 +506 126 38 174 120 +507 223 237 379 388 +508 171 175 379 388 +509 359 368 362 370 +510 211 204 196 370 +511 168 174 38 120 +512 148 156 163 372 +513 152 143 191 385 +514 58 221 165 214 +515 195 189 199 384 +516 58 221 169 165 +517 118 129 189 384 +518 170 371 171 379 +519 191 200 46 152 +520 214 223 237 379 +521 165 175 59 214 +522 46 191 152 143 +523 149 102 94 374 +524 146 363 380 386 +525 206 370 204 375 +526 154 47 191 201 +527 368 370 360 375 +528 104 95 109 375 +529 117 216 35 127 +530 205 206 197 375 +531 362 370 368 375 +532 27 95 190 105 +533 19 107 226 215 +534 54 198 192 166 +535 192 198 204 376 +536 192 203 55 166 +537 19 107 215 96 +538 119 130 140 380 +539 116 100 109 368 +540 94 106 142 374 +541 226 228 215 377 +542 128 145 380 386 +543 101 108 96 377 +544 182 180 176 376 +545 132 120 126 378 +546 101 96 18 224 +547 173 180 167 381 +548 133 369 134 384 +549 220 235 228 373 +550 228 230 220 373 +551 124 140 122 369 +552 216 227 35 127 +553 96 215 18 224 +554 144 156 151 382 +555 187 172 170 371 +556 140 369 133 380 +557 200 152 385 387 +558 205 200 385 387 +559 152 157 385 387 +560 118 42 193 189 +561 181 174 168 378 +562 171 172 188 371 +563 202 105 190 375 +564 178 166 203 376 +565 170 165 169 379 +566 98 100 116 368 +567 198 54 176 166 +568 218 213 217 382 +569 171 175 165 379 +570 218 235 220 373 +571 93 98 97 381 +572 172 171 170 371 +573 128 380 133 386 +574 99 93 103 381 +575 94 14 149 142 +576 146 380 145 386 +577 118 43 189 129 +578 127 123 117 383 +579 122 121 117 383 +580 213 219 225 382 +581 199 189 43 129 +582 237 227 229 388 +583 26 95 104 190 +584 193 125 42 118 +585 191 201 212 385 +586 143 164 154 385 +587 124 133 140 369 +588 133 134 125 384 +589 109 100 110 368 +590 229 227 216 383 +591 193 189 194 384 +592 142 94 15 106 +593 206 196 204 370 +594 205 196 206 370 +595 157 158 148 372 +596 156 148 158 372 +597 203 192 211 376 +598 197 26 104 190 +599 129 134 132 384 +600 102 14 149 94 +601 102 161 113 387 +602 160 207 154 386 +603 202 27 190 105 +604 203 178 55 166 +605 105 116 95 375 +606 131 120 139 378 +607 188 168 179 378 +608 131 138 186 388 +609 107 96 115 377 +610 175 171 179 388 +611 106 110 108 374 +612 156 158 153 374 +613 222 236 216 383 +614 142 106 153 374 +615 155 163 144 382 +616 119 133 128 380 +617 102 149 161 387 +618 154 146 145 386 +619 215 235 224 377 +620 207 201 154 386 +621 179 131 186 388 +622 113 209 104 387 +623 187 177 167 381 +624 204 202 206 375 +625 209 197 104 387 +626 160 48 154 207 +627 161 113 13 102 +628 146 141 145 380 +629 191 205 200 385 +630 143 152 157 385 +631 178 181 182 376 +632 133 124 134 369 +633 200 205 197 387 +634 157 152 149 387 +635 195 194 189 384 +636 186 131 40 138 +637 146 147 141 380 +638 221 214 230 379 +639 136 193 207 386 +640 125 193 136 386 +641 204 211 192 376 +642 165 170 171 379 +643 127 234 138 388 +644 125 128 133 386 +645 95 116 109 375 +646 99 98 93 381 +647 179 131 40 186 +648 232 217 49 151 +649 201 154 48 207 +650 149 161 13 102 +651 21 111 173 97 +652 117 123 122 383 +653 127 227 234 388 +654 145 128 29 160 +655 213 218 219 382 +656 106 15 142 153 +657 159 232 49 151 +658 21 111 185 173 +659 134 129 118 384 +660 160 128 29 136 +661 160 145 128 386 +662 96 108 115 377 +663 139 120 132 378 +664 215 228 235 377 +665 156 144 163 382 +666 160 128 136 386 +667 45 152 200 387 +668 174 184 37 126 +669 199 44 210 135 +670 180 187 167 381 +671 188 181 168 378 +672 94 110 106 374 +673 142 153 158 374 +674 20 114 233 226 +675 104 102 113 387 +676 127 234 36 138 +677 114 20 107 226 +678 216 236 229 383 +679 102 110 94 374 +680 149 142 158 374 +681 185 198 53 176 +682 53 208 198 185 +683 227 234 36 127 +684 17 101 224 232 +685 126 184 37 135 +686 33 231 222 121 +687 17 101 232 112 +688 113 209 25 104 +689 125 41 136 193 +690 41 207 136 193 +691 231 33 137 121 +692 130 150 32 137 +693 137 150 32 162 +694 24 103 177 183 +695 24 103 183 114 +696 202 190 206 375 +697 104 209 25 197 +698 225 52 162 231 +699 112 159 16 153 +700 52 162 155 225 +701 166 178 182 376 +702 169 183 233 57 +703 223 227 237 388 +704 233 221 169 57 +705 140 133 119 380 +706 131 127 138 388 +707 45 200 209 387 +708 152 45 161 387 +709 206 190 197 375 +710 234 223 60 388 +711 208 28 202 111 +712 166 182 176 376 +713 191 212 205 385 +714 164 143 157 385 +715 184 56 210 178 +716 134 118 125 384 +717 105 202 28 111 +718 203 210 56 178 +719 207 193 201 386 +720 44 199 129 135 +721 209 161 45 387 +722 237 230 214 379 +723 13 14 149 102 +724 234 60 186 388 +725 48 47 154 201 +726 217 50 49 151 +727 40 39 131 179 +728 173 22 21 97 +729 51 144 50 213 +730 54 53 198 176 +731 227 36 35 127 +732 18 17 101 224 +733 160 136 207 386 +734 106 112 16 153 +735 93 22 23 167 +736 119 30 31 141 +737 128 29 30 145 +738 125 42 41 193 +739 32 31 150 130 +740 23 24 103 177 +741 34 33 222 121 +742 57 221 169 58 +743 107 20 19 226 +744 104 25 26 197 +745 45 46 200 152 +746 209 113 161 387 +747 165 59 58 214 +748 186 138 234 388 +749 51 52 155 225 +750 37 38 174 126 +751 46 191 143 47 +752 18 19 215 96 +753 168 38 39 120 +754 59 175 60 223 +755 216 35 34 117 +756 207 160 5 136 +757 1 232 159 112 +758 2 209 113 161 +759 111 208 4 185 +760 6 137 162 231 +761 114 183 3 233 +762 184 210 8 135 +763 234 186 7 138 +764 202 28 27 105 +765 199 43 44 129 +766 178 56 203 55 +767 136 128 125 386 +768 145 160 154 386 +769 118 43 42 189 +770 209 200 197 387 +771 152 161 149 387 +772 186 175 179 388 +773 111 4 21 185 +774 1 232 49 159 +775 2 45 209 161 +776 160 29 5 136 +777 60 186 7 234 +778 55 54 192 166 +779 95 26 27 190 +780 159 16 1 112 +781 44 8 210 135 +782 7 36 234 138 +783 52 6 162 231 +784 56 210 8 184 +785 3 20 114 233 +786 186 40 7 138 +787 106 16 15 153 +788 208 53 4 185 +789 207 5 41 136 +790 113 13 2 161 +791 25 2 209 113 +792 137 6 33 231 +793 184 8 37 135 +794 1 17 232 112 +795 48 160 5 207 +796 208 4 28 111 +797 32 6 137 162 +798 183 3 24 114 +799 234 227 223 388 +800 233 183 3 57 +801 15 14 94 142 +802 111 202 362 105 +803 111 362 202 208 +804 178 210 366 203 +805 178 366 210 184 +806 378 124 361 369 +807 378 361 124 139 +808 123 124 361 139 +809 379 228 373 230 +810 379 373 228 364 +811 374 109 368 360 +812 374 368 109 110 +813 123 388 383 127 +814 383 388 123 361 +815 166 198 376 176 +816 166 376 198 192 +817 388 60 175 223 +818 388 175 60 186 +819 122 124 383 123 +820 122 383 124 369 +821 361 383 124 123 +822 361 124 383 369 +3 2 4 701 +823 345 354 296 397 +824 315 257 242 395 +825 251 283 258 393 +826 329 277 281 395 +827 356 346 284 396 +828 273 353 283 393 +829 258 283 353 393 +830 258 353 339 393 +831 257 329 281 395 +832 305 296 354 397 +833 137 356 284 396 +834 257 315 329 395 +835 305 354 138 397 +836 356 137 121 396 +837 137 284 272 396 +838 127 138 354 397 +839 346 275 284 396 +840 329 324 277 395 +841 353 273 347 393 +842 121 343 356 396 +843 339 246 258 393 +844 281 271 257 395 +845 137 272 130 396 +846 281 277 271 395 +847 275 272 284 396 +848 346 356 343 396 +849 249 257 271 395 +850 138 127 131 397 +851 353 347 339 393 +852 258 246 251 393 +853 329 315 324 395 +854 242 257 249 395 +855 276 283 251 393 +856 137 130 121 396 +857 393 399 395 406 +858 395 398 393 406 +859 127 354 348 397 +860 393 398 253 406 +861 389 395 393 399 +862 393 395 389 398 +863 273 283 276 393 +864 345 348 354 397 +865 305 299 296 397 +866 345 397 296 413 +867 315 395 242 412 +868 400 403 392 407 +869 395 399 279 406 +870 260 253 398 406 +871 397 403 400 407 +872 303 392 403 407 +873 139 397 400 407 +874 398 402 393 411 +875 391 398 254 411 +876 275 396 346 409 +877 400 402 397 414 +878 396 402 400 414 +879 391 402 398 411 +880 351 402 396 414 +881 255 254 245 411 +882 392 400 132 404 +883 303 403 392 415 +884 139 397 123 400 +885 326 401 392 404 +886 277 395 324 408 +887 389 392 400 403 +888 278 399 393 406 +889 347 393 273 409 +890 397 400 389 402 +891 396 402 351 409 +892 279 399 395 408 +893 254 398 245 411 +894 132 400 392 407 +895 310 403 303 415 +896 393 398 389 402 +897 394 405 244 416 +898 396 389 400 402 +899 272 130 396 410 +900 339 246 393 411 +901 322 327 390 408 +902 139 400 132 407 +903 251 393 253 406 +904 244 412 394 416 +905 392 401 326 415 +906 133 400 390 404 +907 390 401 327 404 +908 389 400 397 403 +909 260 253 245 398 +910 260 398 395 406 +911 279 395 277 408 +912 394 398 244 405 +913 346 396 351 409 +914 244 398 394 412 +915 392 401 400 404 +916 327 322 390 419 +917 279 271 277 395 +918 350 397 402 414 +919 333 326 401 415 +920 346 343 351 396 +921 121 343 396 414 +922 322 408 390 419 +923 252 300 405 417 +924 391 389 398 402 +925 391 405 300 417 +926 278 393 399 409 +927 286 390 399 408 +928 392 394 389 401 +929 393 402 399 409 +930 251 246 253 393 +931 389 394 392 403 +932 399 402 396 409 +933 243 395 260 398 +934 333 401 394 415 +935 310 394 403 415 +936 393 358 402 409 +937 400 401 390 404 +938 285 399 396 409 +939 280 279 399 406 +940 132 392 126 407 +941 350 402 397 413 +942 293 300 391 405 +943 357 391 402 413 +944 276 393 251 406 +945 285 396 275 409 +946 310 394 292 403 +947 317 394 333 401 +948 249 395 271 406 +949 273 393 278 409 +950 126 407 392 418 +951 344 391 357 413 +952 254 255 262 411 +953 334 395 401 408 +954 268 286 390 399 +955 327 319 390 404 +956 133 390 125 404 +957 303 403 304 407 +958 254 398 391 405 +959 139 131 397 407 +960 140 400 396 410 +961 310 294 303 403 +962 350 397 345 413 +963 296 397 302 413 +964 309 391 297 413 +965 340 393 358 402 +966 389 400 392 401 +967 274 390 286 408 +968 285 269 396 399 +969 357 391 341 402 +970 324 395 334 408 +971 339 393 340 411 +972 242 395 243 412 +973 316 395 315 412 +974 358 351 402 409 +975 249 260 395 406 +976 254 262 391 411 +977 391 309 403 413 +978 132 126 392 418 +979 341 357 344 391 +980 286 399 279 408 +981 347 358 393 409 +982 275 269 272 396 +983 285 269 275 396 +984 123 400 397 414 +985 351 352 402 414 +986 140 396 130 410 +987 271 395 279 406 +988 241 252 300 405 +989 127 123 131 397 +990 334 395 316 401 +991 125 390 133 419 +992 303 295 392 407 +993 122 396 140 400 +994 269 272 396 410 +995 131 123 139 397 +996 309 293 297 391 +997 396 399 389 402 +998 347 340 339 393 +999 393 389 399 402 +1000 242 249 243 395 +1001 316 324 315 395 +1002 300 287 241 405 +1003 309 293 391 403 +1004 311 315 242 412 +1005 253 393 246 411 +1006 268 286 274 390 +1007 293 391 300 417 +1008 121 130 122 396 +1009 303 392 295 418 +1010 243 249 260 395 +1011 334 324 316 395 +1012 238 311 242 412 +1013 358 340 347 393 +1014 302 403 397 407 +1015 139 124 132 400 +1016 121 396 122 414 +1017 318 333 326 401 +1018 140 122 130 396 +1019 127 397 348 414 +1020 322 274 408 419 +1021 317 394 401 412 +1022 327 390 319 419 +1023 320 301 415 416 +1024 292 403 394 405 +1025 125 404 390 419 +1026 292 291 405 416 +1027 127 123 397 414 +1028 341 391 344 417 +1029 327 401 390 408 +1030 134 132 400 404 +1031 326 328 401 404 +1032 133 390 400 410 +1033 278 399 285 409 +1034 133 134 400 404 +1035 254 391 252 405 +1036 328 327 401 404 +1037 351 396 343 414 +1038 119 267 128 410 +1039 345 296 338 413 +1040 398 403 391 405 +1041 391 413 344 417 +1042 392 407 295 418 +1043 131 299 397 407 +1044 278 280 399 406 +1045 276 278 273 393 +1046 261 247 244 405 +1047 349 248 239 411 +1048 296 290 338 413 +1049 389 391 398 403 +1050 350 357 402 413 +1051 394 415 301 416 +1052 252 405 391 417 +1053 350 348 345 397 +1054 302 296 299 397 +1055 320 415 394 416 +1056 255 253 246 411 +1057 297 391 293 417 +1058 394 403 398 405 +1059 268 390 274 419 +1060 126 295 407 418 +1061 325 238 250 412 +1062 247 244 405 416 +1063 297 293 300 417 +1064 399 401 395 408 +1065 267 119 263 410 +1066 389 390 400 401 +1067 301 289 314 415 +1068 292 405 394 416 +1069 311 238 325 412 +1070 123 122 400 414 +1071 342 358 351 402 +1072 252 391 254 417 +1073 129 392 132 404 +1074 326 392 323 404 +1075 316 401 395 412 +1076 390 404 319 419 +1077 286 270 279 399 +1078 395 399 389 401 +1079 350 402 352 414 +1080 390 408 274 419 +1081 394 401 392 415 +1082 391 403 293 405 +1083 392 403 394 415 +1084 389 398 394 403 +1085 338 290 297 413 +1086 327 319 322 419 +1087 268 399 390 410 +1088 74 311 315 242 +1089 248 411 349 417 +1090 122 396 400 414 +1091 323 404 392 418 +1092 270 278 280 399 +1093 267 128 410 419 +1094 350 342 352 402 +1095 390 401 399 408 +1096 292 293 403 405 +1097 139 123 124 400 +1098 278 393 276 406 +1099 297 413 391 417 +1100 317 333 318 401 +1101 74 311 242 238 +1102 316 317 401 412 +1103 323 392 326 418 +1104 310 292 294 403 +1105 261 244 398 405 +1106 397 403 302 413 +1107 397 402 389 403 +1108 297 344 413 417 +1109 333 394 320 415 +1110 340 393 402 411 +1111 269 396 399 410 +1112 402 403 397 413 +1113 247 287 291 405 +1114 299 302 397 407 +1115 341 402 391 411 +1116 337 275 346 409 +1117 395 398 243 412 +1118 239 335 349 411 +1119 254 391 262 417 +1120 272 263 130 410 +1121 322 274 313 408 +1122 320 333 317 394 +1123 390 410 133 419 +1124 266 347 273 409 +1125 119 130 263 410 +1126 277 324 265 408 +1127 287 247 241 405 +1128 338 290 91 297 +1129 339 335 246 411 +1130 342 351 352 402 +1131 314 320 301 415 +1132 250 325 75 238 +1133 132 392 129 418 +1134 247 405 291 416 +1135 254 261 398 405 +1136 325 412 250 416 +1137 268 410 390 419 +1138 390 389 399 401 +1139 270 280 279 399 +1140 126 295 120 407 +1141 301 394 310 415 +1142 297 344 338 413 +1143 292 394 301 416 +1144 295 288 120 407 +1145 121 117 336 414 +1146 278 270 285 399 +1147 311 75 325 238 +1148 131 120 288 407 +1149 125 319 404 419 +1150 348 397 350 414 +1151 335 239 246 411 +1152 134 132 124 400 +1153 312 118 125 404 +1154 292 310 301 394 +1155 303 294 304 403 +1156 391 403 402 413 +1157 389 394 398 401 +1158 289 87 301 314 +1159 326 328 318 401 +1160 350 357 342 402 +1161 274 265 313 408 +1162 291 292 301 416 +1163 389 398 395 401 +1164 129 392 404 418 +1165 337 266 275 409 +1166 128 119 30 267 +1167 327 334 401 408 +1168 265 324 313 408 +1169 287 70 247 291 +1170 391 389 402 403 +1171 296 345 338 90 +1172 128 133 410 419 +1173 317 320 394 416 +1174 398 401 394 412 +1175 347 266 337 409 +1176 133 400 140 410 +1177 70 241 287 247 +1178 341 411 391 417 +1179 303 415 392 418 +1180 296 338 290 90 +1181 124 123 122 400 +1182 344 338 91 297 +1183 293 287 300 405 +1184 268 267 410 419 +1185 121 336 343 414 +1186 389 390 399 400 +1187 267 119 30 263 +1188 302 304 403 407 +1189 130 31 119 263 +1190 134 133 125 404 +1191 327 328 319 404 +1192 392 415 326 418 +1193 391 411 262 417 +1194 316 317 318 401 +1195 395 401 398 412 +1196 294 292 293 403 +1197 31 130 272 263 +1198 389 399 396 400 +1199 249 264 240 406 +1200 394 412 317 416 +1201 269 399 268 410 +1202 248 262 411 417 +1203 132 129 126 418 +1204 321 314 298 415 +1205 341 349 411 417 +1206 87 320 301 314 +1207 339 246 335 66 +1208 309 302 403 413 +1209 246 335 66 239 +1210 249 271 264 406 +1211 298 314 289 415 +1212 340 402 341 411 +1213 248 262 239 411 +1214 86 289 298 314 +1215 321 86 298 314 +1216 125 319 312 404 +1217 240 251 253 406 +1218 241 254 252 405 +1219 325 317 412 416 +1220 250 412 244 416 +1221 399 400 390 410 +1222 243 260 245 398 +1223 280 271 279 406 +1224 120 39 131 288 +1225 245 254 261 398 +1226 303 304 295 407 +1227 346 337 82 275 +1228 67 349 248 239 +1229 396 400 399 410 +1230 132 126 120 407 +1231 322 313 327 408 +1232 240 276 251 406 +1233 277 265 324 79 +1234 274 322 313 78 +1235 279 277 265 408 +1236 323 118 312 404 +1237 337 82 275 266 +1238 337 346 351 409 +1239 266 273 278 409 +1240 119 128 133 410 +1241 133 124 140 400 +1242 238 242 243 412 +1243 311 316 315 412 +1244 349 67 335 239 +1245 324 265 313 79 +1246 313 274 78 265 +1247 273 266 83 347 +1248 327 334 318 401 +1249 127 348 117 414 +1250 83 266 337 347 +1251 350 345 338 413 +1252 296 302 290 413 +1253 117 348 336 414 +1254 268 270 286 399 +1255 340 358 342 402 +1256 336 34 117 121 +1257 240 264 276 406 +1258 285 270 269 399 +1259 339 340 335 411 +1260 71 241 300 287 +1261 351 343 352 414 +1262 131 288 299 407 +1263 244 243 398 412 +1264 321 415 298 418 +1265 321 326 415 418 +1266 252 254 262 417 +1267 342 357 341 402 +1268 269 270 268 399 +1269 263 268 267 410 +1270 122 140 124 400 +1271 309 294 293 403 +1272 314 333 320 415 +1273 269 263 272 410 +1274 129 404 323 418 +1275 334 316 318 401 +1276 302 304 294 403 +1277 127 117 123 414 +1278 118 125 42 312 +1279 298 415 303 418 +1280 341 349 335 411 +1281 121 122 117 414 +1282 340 342 341 402 +1283 294 309 302 403 +1284 121 336 34 343 +1285 300 252 241 71 +1286 250 238 244 412 +1287 325 317 311 412 +1288 292 291 287 405 +1289 301 310 289 415 +1290 133 134 124 400 +1291 328 327 318 401 +1292 62 264 240 249 +1293 323 118 43 312 +1294 244 243 245 398 +1295 336 117 35 348 +1296 348 117 35 127 +1297 299 131 39 288 +1298 139 120 131 407 +1299 312 125 42 319 +1300 261 244 245 398 +1301 249 62 264 271 +1302 276 63 240 251 +1303 331 125 136 419 +1304 249 240 260 406 +1305 344 357 338 413 +1306 290 309 297 413 +1307 240 63 276 264 +1308 134 129 132 404 +1309 326 323 328 404 +1310 285 275 266 409 +1311 129 118 323 404 +1312 320 332 306 416 +1313 297 308 355 417 +1314 267 268 274 419 +1315 274 286 265 408 +1316 347 337 358 409 +1317 313 324 334 408 +1318 331 319 125 419 +1319 241 247 261 405 +1320 119 140 130 410 +1321 297 355 344 417 +1322 331 282 322 419 +1323 301 320 306 416 +1324 250 244 247 416 +1325 317 325 320 416 +1326 322 282 274 419 +1327 349 341 344 417 +1328 125 136 41 331 +1329 61 281 271 257 +1330 278 276 280 406 +1331 299 304 302 407 +1332 122 123 117 414 +1333 256 291 306 416 +1334 256 247 291 416 +1335 246 239 255 411 +1336 126 135 307 418 +1337 355 92 308 297 +1338 295 126 307 418 +1339 298 289 303 415 +1340 259 308 300 417 +1341 139 132 120 407 +1342 321 326 314 415 +1343 292 287 293 405 +1344 332 320 306 88 +1345 269 268 263 410 +1346 128 125 133 419 +1347 240 253 260 406 +1348 89 345 354 296 +1349 298 303 295 418 +1350 315 242 257 73 +1351 259 300 252 417 +1352 352 348 350 414 +1353 261 254 241 405 +1354 92 355 344 297 +1355 319 125 41 331 +1356 271 249 61 257 +1357 340 341 335 411 +1358 250 332 76 325 +1359 301 306 320 88 +1360 316 311 317 412 +1361 244 238 243 412 +1362 118 323 43 129 +1363 305 89 354 296 +1364 315 257 329 73 +1365 266 278 285 409 +1366 321 323 326 418 +1367 330 135 129 418 +1368 334 327 313 408 +1369 357 350 338 413 +1370 290 302 309 413 +1371 279 265 286 408 +1372 332 325 250 416 +1373 351 358 337 409 +1374 256 332 76 250 +1375 251 258 283 64 +1376 85 321 298 418 +1377 129 134 118 404 +1378 312 328 323 404 +1379 119 133 140 410 +1380 128 267 29 419 +1381 69 247 291 256 +1382 138 36 127 354 +1383 306 69 291 256 +1384 356 346 81 284 +1385 284 346 81 275 +1386 319 331 322 419 +1387 250 256 332 416 +1388 300 308 297 417 +1389 125 118 134 404 +1390 328 312 319 404 +1391 295 126 37 307 +1392 347 84 273 353 +1393 307 126 37 135 +1394 84 283 273 353 +1395 252 262 248 417 +1396 33 121 343 356 +1397 121 33 137 356 +1398 281 277 329 80 +1399 322 331 77 282 +1400 259 308 72 300 +1401 137 32 130 272 +1402 32 137 284 272 +1403 324 329 277 80 +1404 77 274 322 282 +1405 44 129 330 135 +1406 258 353 65 339 +1407 280 276 264 406 +1408 339 65 258 246 +1409 288 304 299 407 +1410 300 72 259 252 +1411 298 307 85 418 +1412 264 271 280 406 +1413 305 138 40 299 +1414 136 128 29 419 +1415 330 321 85 418 +1416 295 304 288 407 +1417 282 29 267 419 +1418 131 40 138 299 +1419 248 68 259 417 +1420 336 348 352 414 +1421 326 333 314 415 +1422 289 310 303 415 +1423 291 301 306 416 +1424 330 129 323 418 +1425 352 343 336 414 +1426 276 251 283 64 +1427 311 74 75 238 +1428 307 330 85 418 +1429 136 29 282 419 +1430 91 92 344 297 +1431 255 239 262 411 +1432 36 348 127 354 +1433 41 42 125 319 +1434 296 89 345 90 +1435 271 249 62 61 +1436 88 320 301 87 +1437 242 315 74 73 +1438 129 135 126 418 +1439 355 259 68 417 +1440 70 247 291 69 +1441 86 298 85 321 +1442 126 37 38 295 +1443 83 84 273 347 +1444 330 129 44 323 +1445 31 30 119 263 +1446 30 29 128 267 +1447 256 306 332 416 +1448 121 34 33 343 +1449 130 32 31 272 +1450 246 339 65 66 +1451 346 82 81 275 +1452 79 324 277 80 +1453 78 274 322 77 +1454 282 331 136 419 +1455 307 135 330 418 +1456 67 66 335 239 +1457 308 259 355 417 +1458 252 300 72 71 +1459 90 338 290 91 +1460 289 86 87 314 +1461 337 82 266 83 +1462 265 78 313 79 +1463 287 71 241 70 +1464 10 281 257 329 +1465 259 355 11 308 +1466 330 307 8 135 +1467 9 258 283 353 +1468 6 137 356 284 +1469 256 306 12 332 +1470 331 136 5 282 +1471 354 138 7 305 +1472 76 250 325 75 +1473 248 349 67 68 +1474 39 40 131 299 +1475 320 325 332 416 +1476 250 247 256 416 +1477 251 276 63 64 +1478 128 136 125 419 +1479 307 298 295 418 +1480 344 355 349 417 +1481 282 267 274 419 +1482 35 36 348 127 +1483 136 29 5 282 +1484 68 355 11 259 +1485 329 257 10 73 +1486 89 354 7 305 +1487 85 307 8 330 +1488 12 256 332 76 +1489 240 63 264 62 +1490 117 35 34 336 +1491 38 39 120 288 +1492 283 258 9 64 +1493 8 44 330 135 +1494 11 259 308 72 +1495 138 7 36 354 +1496 40 7 138 305 +1497 81 6 356 284 +1498 306 12 69 256 +1499 355 92 11 308 +1500 129 43 44 323 +1501 248 259 252 417 +1502 330 323 321 418 +1503 33 6 137 356 +1504 5 41 136 331 +1505 307 37 8 135 +1506 329 10 281 80 +1507 10 281 61 257 +1508 84 9 283 353 +1509 332 306 12 88 +1510 77 331 5 282 +1511 137 6 32 284 +1512 353 9 258 65 +1513 118 42 43 312 +1514 393 411 245 398 +1515 245 253 393 398 +1516 411 253 245 255 +1517 411 245 253 393 +1518 299 138 397 305 +1519 299 397 138 131 +1520 120 295 38 288 +1521 38 295 120 126 +1522 417 68 349 248 +1523 417 349 68 355 +$EndElements diff --git a/docs/examples/meshes/interface.msh b/docs/examples/meshes/interface.msh new file mode 100644 index 000000000..deb7fa2d1 --- /dev/null +++ b/docs/examples/meshes/interface.msh @@ -0,0 +1,432 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +2 +1 9 "interfacee" +2 12 "both" +$EndPhysicalNames +$Entities +6 7 2 0 +1 -0.5 0 0 0 +2 0.5 0 0 0 +3 0.5 0.5 0 0 +4 -0.5 0.5 0 0 +5 0.5 1 0 0 +6 -0.5 1 0 0 +1 -0.5000000999999999 -1e-07 -1e-07 0.5000000999999999 1e-07 1e-07 0 2 1 -2 +2 0.4999999 -9.999999997511999e-08 -1e-07 0.5000000999999999 0.5000000999999999 1e-07 0 2 2 -3 +3 -0.5000001000000001 0.4999998999999999 -1.000000000555111e-07 0.5000001000000001 0.5000001000000001 1.000000000555111e-07 1 9 2 3 -4 +4 -0.5000000999999999 -9.999999997511999e-08 -1e-07 -0.4999999 0.5000000999999999 1e-07 0 2 4 -1 +5 0.4999999 0.4999999 -1e-07 0.5000000999999999 1.0000001 1e-07 0 2 3 -5 +6 -0.5000000999999999 0.9999999000000001 -1e-07 0.5000000999999999 1.0000001 1e-07 0 2 5 -6 +7 -0.5000000999999999 0.4999999 -1e-07 -0.4999999 1.0000001 1e-07 0 2 6 -4 +1 -0.5000001000000001 -1.000000000583867e-07 -1.000000000555111e-07 0.5000001000000001 0.5000001000000001 1.000000000555111e-07 1 12 4 1 2 3 4 +2 -0.5000001000000001 0.4999998999999999 -1.000000000555111e-07 0.5000001000000001 1.0000001 1.000000000555111e-07 1 12 4 3 -7 -6 -5 +$EndEntities +$Nodes +15 102 1 102 +0 1 0 1 +1 +-0.5 0 0 +0 2 0 1 +2 +0.5 0 0 +0 3 0 1 +3 +0.5 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.5 1 0 +1 1 0 7 +7 +8 +9 +10 +11 +12 +13 +-0.3750000000000002 0 0 +-0.2500000000000002 0 0 +-0.1250000000000006 0 0 +-1.165734175856414e-15 0 0 +0.1249999999999988 0 0 +0.2499999999999991 0 0 +0.3749999999999994 0 0 +1 2 0 3 +14 +15 +16 +0.5 0.1249999999999999 0 +0.5 0.2499999999999994 0 +0.5 0.3749999999999996 0 +1 3 0 7 +17 +18 +19 +20 +21 +22 +23 +0.3750000000000002 0.5 0 +0.2500000000000002 0.5 0 +0.1250000000000006 0.5 0 +1.165734175856414e-15 0.5 0 +-0.1249999999999988 0.5 0 +-0.2499999999999991 0.5 0 +-0.3749999999999994 0.5 0 +1 4 0 3 +24 +25 +26 +-0.5 0.3750000000000001 0 +-0.5 0.2500000000000006 0 +-0.5 0.1250000000000004 0 +1 5 0 3 +27 +28 +29 +0.5 0.6249999999999999 0 +0.5 0.7499999999999994 0 +0.5 0.8749999999999996 0 +1 6 0 7 +30 +31 +32 +33 +34 +35 +36 +0.3750000000000002 1 0 +0.2500000000000002 1 0 +0.1250000000000006 1 0 +1.165734175856414e-15 1 0 +-0.1249999999999988 1 0 +-0.2499999999999991 1 0 +-0.3749999999999994 1 0 +1 7 0 3 +37 +38 +39 +-0.5 0.8750000000000001 0 +-0.5 0.7500000000000006 0 +-0.5 0.6250000000000004 0 +2 1 0 32 +40 +41 +42 +43 +44 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +-0.06250000000000089 0.1082531754730544 0 +0.1886995265858526 0.3869838548559241 0 +-0.1864783727984019 0.3923366612667986 0 +0.1864783727984019 0.1076633387332014 0 +-0.40127596735477 0.1791127824134674 0 +0.3800901487253206 0.1987663973443807 0 +-0.06232972879973279 0.3918451306502542 0 +-0.1248013502663551 0.288904729145801 0 +-0.0103551798443477 0.2965859159320229 0 +-0.2485292786755291 0.2934157229184406 0 +-0.1874999999999988 0.1752404735808357 0 +-0.3856768717468466 0.3059427571357878 0 +0.06800036305122908 0.1048809271381685 0 +0.1361631226416047 0.2036719675913905 0 +0.2495138652037691 0.2008441225213108 0 +0.3196600551454468 0.3998525713990545 0 +0.05533180207213279 0.3964253672935926 0 +0.4014541938652549 0.3215633842485699 0 +0.3077219766759978 0.1012016882888624 0 +-0.3090296953791111 0.4000314650624863 0 +-0.3006692781366349 0.1025606614856545 0 +-0.06644211475527766 0.2101568659081853 0 +0.1931668038378242 0.278380986997215 0 +-0.1875000000000004 0.07647483778475356 0 +-0.2969994706343269 0.2001385014571439 0 +-0.4084936490538903 0.4084936490538906 0 +0.4084936490538905 0.4084936490538904 0 +0.4084936490538905 0.09150635094610934 0 +-0.4114110529975507 0.08739497392595721 0 +0.2921118735719 0.2975992580752718 0 +0.09532761573727341 0.2980426574599756 0 +0.02879177857762211 0.2016740457055651 0 +2 2 0 31 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +86 +87 +88 +89 +90 +91 +92 +93 +94 +95 +96 +97 +98 +99 +100 +101 +102 +0.04950340078699763 0.6104767907801998 0 +-0.1874999999999989 0.6082531754730549 0 +0.1875000000000004 0.8970432141416622 0 +-0.0624999999999988 0.8917468245269453 0 +-0.3119703288398897 0.8941534943680208 0 +0.3727745823942785 0.6922959892733603 0 +-0.4011617204020835 0.6795950867083989 0 +-0.187411721473314 0.8921479361671245 0 +-0.1249852869122179 0.7888568906086373 0 +0.009850869244568227 0.7826545169346127 0 +-0.06249999999999828 0.6752404735808362 0 +-0.250662295632269 0.7980878577805447 0 +0.3119703288398905 0.8944751895584195 0 +0.2524004431642587 0.7867072143630472 0 +0.06414181154076204 0.897786090390902 0 +0.1841909972544986 0.6027808049749023 0 +0.1250000000000006 0.7165063509461087 0 +0.3921815513536758 0.8145548345887761 0 +-0.3198293268430165 0.5965557317679157 0 +0.3066219291087071 0.5952617797139239 0 +-0.3881508355468093 0.8035071324068801 0 +0.127778624789918 0.8161394773552666 0 +-0.06249999999999881 0.5764748377847537 0 +-0.1853525963066585 0.7120651139230396 0 +0.2481975903843487 0.6787104278542685 0 +-0.40849364905389 0.9084936490538902 0 +-0.4084936490538902 0.5915063509461098 0 +0.4101976870872525 0.585081495839929 0 +0.4084936490538904 0.9084936490538904 0 +-0.2947265520143724 0.6999010530413391 0 +0.03402380550897779 0.6968626976766455 0 +$EndNodes +$Elements +3 178 19 216 +1 3 1 8 +19 3 17 +20 17 18 +21 18 19 +22 19 20 +23 20 21 +24 21 22 +25 22 23 +26 23 4 +2 1 2 86 +47 51 64 49 +48 45 69 54 +49 44 64 51 +50 54 58 45 +51 49 59 51 +52 41 69 55 +53 60 63 50 +54 62 69 41 +55 40 61 50 +56 50 63 40 +57 50 64 60 +58 21 46 20 +59 48 56 46 +60 10 52 40 +61 11 52 10 +62 46 56 20 +63 42 46 21 +64 26 44 25 +65 44 51 25 +66 12 43 11 +67 49 50 47 +68 42 47 46 +69 42 49 47 +70 43 52 11 +71 22 42 21 +72 25 51 24 +73 43 54 53 +74 47 48 46 +75 43 53 52 +76 10 40 9 +77 18 55 17 +78 15 45 14 +79 19 41 18 +80 20 56 19 +81 19 56 41 +82 41 55 18 +83 57 69 45 +84 16 57 15 +85 13 58 12 +86 42 59 49 +87 23 59 22 +88 43 58 54 +89 15 57 45 +90 12 58 43 +91 22 59 42 +92 8 60 7 +93 50 61 47 +94 40 63 9 +95 3 66 16 +96 24 65 4 +97 17 66 3 +98 4 65 23 +99 14 67 2 +100 2 67 13 +101 52 71 40 +102 56 70 41 +103 1 68 26 +104 7 68 1 +105 8 63 60 +106 40 71 61 +107 41 70 62 +108 49 64 50 +109 13 67 58 +110 23 65 59 +111 16 66 57 +112 54 62 53 +113 47 61 48 +114 9 63 8 +115 26 68 44 +116 51 65 24 +117 45 67 14 +118 55 66 17 +119 48 70 56 +120 53 71 52 +121 58 67 45 +122 59 65 51 +123 57 66 55 +124 60 64 44 +125 44 68 60 +126 60 68 7 +127 48 71 70 +128 70 71 53 +129 55 69 57 +130 54 69 62 +131 61 71 48 +132 62 70 53 +2 2 2 84 +133 77 89 85 +134 85 96 77 +135 92 101 83 +136 76 92 83 +137 85 89 84 +138 73 94 82 +139 85 93 88 +140 82 95 73 +141 88 93 81 +142 82 94 72 +143 88 96 85 +144 78 101 92 +145 34 79 75 +146 79 83 80 +147 35 79 34 +148 79 80 75 +149 75 86 33 +150 36 76 35 +151 76 79 35 +152 76 83 79 +153 18 87 19 +154 32 74 31 +155 80 81 75 +156 81 86 75 +157 34 75 33 +158 19 87 72 +159 21 73 22 +160 39 78 38 +161 80 82 81 +162 31 84 30 +163 28 77 27 +164 32 86 74 +165 19 72 20 +166 87 88 72 +167 74 85 84 +168 74 84 31 +169 33 86 32 +170 77 96 91 +171 73 101 90 +172 81 102 88 +173 95 101 73 +174 72 102 82 +175 22 90 23 +176 17 91 18 +177 29 89 28 +178 38 92 37 +179 28 89 77 +180 73 90 22 +181 18 91 87 +182 78 92 38 +183 80 95 82 +184 72 94 20 +185 21 94 73 +186 74 93 85 +187 81 93 86 +188 6 97 36 +189 37 97 6 +190 5 100 29 +191 30 100 5 +192 4 98 39 +193 23 98 4 +194 27 99 3 +195 3 99 17 +196 87 96 88 +197 90 98 23 +198 17 99 91 +199 29 100 89 +200 92 97 37 +201 83 95 80 +202 86 93 74 +203 20 94 21 +204 36 97 76 +205 84 100 30 +206 77 99 27 +207 39 98 78 +208 82 102 81 +209 88 102 72 +210 91 99 77 +211 89 100 84 +212 78 98 90 +213 76 97 92 +214 91 96 87 +215 90 101 78 +216 83 101 95 +$EndElements diff --git a/docs/examples/meshes/internal.msh b/docs/examples/meshes/internal.msh new file mode 100644 index 000000000..8166aa7c9 --- /dev/null +++ b/docs/examples/meshes/internal.msh @@ -0,0 +1,686 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +6 +1 7 "top" +1 8 "bottom" +1 9 "left" +1 10 "right" +1 11 "internal" +2 6 "domain" +$EndPhysicalNames +$Entities +6 5 1 0 +1 -0.5 -0.5 0 0 +2 0.5 -0.5 0 0 +3 0.5 0.5 0 0 +4 -0.5 0.5 0 0 +5 0.1 0.1 0 0 +6 0.4 0.4 0 0 +1 -0.5 -0.5 0 0.5 -0.5 0 1 8 2 1 -2 +2 0.5 -0.5 0 0.5 0.5 0 1 10 2 2 -3 +3 -0.5 0.5 0 0.5 0.5 0 1 7 2 3 -4 +4 -0.5 -0.5 0 -0.5 0.5 0 1 9 2 4 -1 +5 0.1 0.1 0 0.4 0.4 0 1 11 2 5 -6 +1 -0.5 -0.5 0 0.5 0.5 0 1 6 4 3 4 1 2 +$EndEntities +$Nodes +12 158 1 158 +0 1 0 1 +1 +-0.5 -0.5 0 +0 2 0 1 +2 +0.5 -0.5 0 +0 3 0 1 +3 +0.5 0.5 0 +0 4 0 1 +4 +-0.5 0.5 0 +0 5 0 1 +5 +0.1 0.1 0 +0 6 0 1 +6 +0.4 0.4 0 +1 1 0 9 +7 +8 +9 +10 +11 +12 +13 +14 +15 +-0.4000000000002774 -0.5 0 +-0.3000000000005548 -0.5 0 +-0.2000000000008322 -0.5 0 +-0.1000000000011096 -0.5 0 +-1.376121439022882e-12 -0.5 0 +0.09999999999889042 -0.5 0 +0.1999999999991677 -0.5 0 +0.2999999999994452 -0.5 0 +0.3999999999997227 -0.5 0 +1 2 0 9 +16 +17 +18 +19 +20 +21 +22 +23 +24 +0.5 -0.4000000000002774 0 +0.5 -0.3000000000005548 0 +0.5 -0.2000000000008322 0 +0.5 -0.1000000000011096 0 +0.5 -1.376121439022882e-12 0 +0.5 0.09999999999889042 0 +0.5 0.1999999999991677 0 +0.5 0.2999999999994452 0 +0.5 0.3999999999997227 0 +1 3 0 9 +25 +26 +27 +28 +29 +30 +31 +32 +33 +0.4000000000002774 0.5 0 +0.3000000000005548 0.5 0 +0.2000000000008322 0.5 0 +0.1000000000011096 0.5 0 +1.376121439022882e-12 0.5 0 +-0.09999999999889042 0.5 0 +-0.1999999999991677 0.5 0 +-0.2999999999994452 0.5 0 +-0.3999999999997227 0.5 0 +1 4 0 9 +34 +35 +36 +37 +38 +39 +40 +41 +42 +-0.5 0.4000000000002774 0 +-0.5 0.3000000000005548 0 +-0.5 0.2000000000008322 0 +-0.5 0.1000000000011096 0 +-0.5 1.376121439022882e-12 0 +-0.5 -0.09999999999889042 0 +-0.5 -0.1999999999991677 0 +-0.5 -0.2999999999994452 0 +-0.5 -0.3999999999997227 0 +1 5 0 4 +43 +44 +45 +46 +0.15999999999987 0.15999999999987 0 +0.2199999999996804 0.2199999999996804 0 +0.279999999999662 0.279999999999662 0 +0.33999999999983 0.33999999999983 0 +2 1 0 112 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +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 +-0.4243412329254226 -0.1631229621210719 0 +-0.04890594039226279 -0.413609702691484 0 +-0.4164293427566662 0.142764687502224 0 +0.4037434561857105 -0.1547726937312097 0 +-0.1392918823219163 0.4091564794319314 0 +0.1493803741801924 -0.4136716522835989 0 +-0.2469062977956384 -0.4129247352618338 0 +0.41493091456133 0.05249669290835797 0 +0.04737409167734917 0.4095676866458485 0 +-0.3498327514899938 0.421672967468681 0 +0.3453499396116337 -0.4147990284844523 0 +-0.4073668444074327 -0.3402905459709372 0 +0.4202480801529463 0.2413301707333462 0 +0.2439069846004146 0.4252733639538581 0 +0.4163799523371253 0.1488575493878978 0 +0.3329474309357859 0.1010408251074908 0 +0.3305680314384424 0.007504560189976516 0 +0.2521921749178729 0.05460050440420325 0 +0.2476644651500318 -0.0345803547071032 0 +0.1595796595911618 0.02030602955348995 0 +0.1651351597447345 -0.0751434341441332 0 +0.08919633395596191 -0.02852573276507249 0 +0.08821750500287141 -0.1151762244872123 0 +0.01014397052338416 -0.06971258826452294 0 +0.005285398200761075 -0.1658333030593626 0 +-0.07369878487653934 -0.1092408281507412 0 +-0.06760869567442712 -0.01931222833559496 0 +-0.1491656956435887 -0.06049652249657559 0 +-0.144887539646979 0.03102347157211024 0 +-0.1497814059214502 -0.1449210091650885 0 +0.1636103584871826 -0.1668651826375512 0 +-0.2272504652766112 -0.01066731738100382 0 +-0.2234453449591219 0.08168451222636204 0 +-0.139636507870723 0.124309783955579 0 +-0.2189333177904932 0.1757499833959893 0 +-0.134392601640576 0.2182809294609523 0 +-0.05716275847972697 0.1665533746008246 0 +-0.05060700848188666 0.2601163616706027 0 +-0.2130442590886793 0.2679002671491382 0 +-0.3063131594077195 0.03611441523939263 0 +0.02362507177378339 0.2107272377986898 0 +-0.2958517905218295 0.2297544262758539 0 +-0.3117620292332341 -0.05665176693343683 0 +0.03517792816966737 0.3049319077271943 0 +-0.04480587224024518 0.3534124380787279 0 +0.1245846178001701 0.3438456353738175 0 +0.323997935550016 -0.08437925933011475 0 +0.3267850298182902 0.199679572268237 0 +-0.3941470033768285 -4.030463711477078e-05 0 +0.1214929250220547 0.2446345604981789 0 +-0.3017394530211586 0.3316758928311875 0 +-0.3820063117624554 0.2796381080466135 0 +0.0501551264921184 -0.4140143110609256 0 +0.09944040600223227 -0.3308168546400662 0 +0.1989036021284097 -0.3301096860855161 0 +0.0009667559064284829 -0.3295861673299127 0 +-0.09630784331434823 -0.3285798504170381 0 +-0.04453323654267555 -0.2499718366837912 0 +-0.1476614291693811 -0.4128962931366361 0 +-0.1946088958612251 -0.3259412844556226 0 +-0.2962532787961263 -0.328879143836835 0 +-0.2472503589364749 -0.2444207815985844 0 +-0.3490026389935644 -0.2282949502243351 0 +-0.2845138377984643 -0.1527374930500327 0 +0.2489324922693545 -0.4142252919745167 0 +0.2980625370521982 -0.3273862427456999 0 +0.2443583319432001 -0.2297956380336114 0 +0.3464627853628836 -0.2407883903108089 0 +0.4065444220978809 -0.3364235125566677 0 +-0.2485848778447144 0.418730672382845 0 +0.4146791428034675 -0.04778005961083855 0 +-0.06347101078338376 0.07316565584849184 0 +0.05083473522943732 -0.253383882324431 0 +0.01019585781482664 0.01792716872978516 0 +-0.3149650670534511 0.1343319192025499 0 +0.1980384757726949 0.3019615242266475 0 +-0.1295921237804424 0.3124046545542715 0 +0.2427241987162503 -0.123916557531307 0 +0.02578903471310925 0.110958264417823 0 +-0.1487551827215856 -0.235708991856202 0 +0.2508402914347567 0.1399660790869116 0 +-0.3462202233738321 -0.4181482240520701 0 +0.1453559389007542 -0.2521463033751026 0 +0.1530287302008692 0.4244150892092525 0 +-0.2268348248073812 -0.09736195031593869 0 +-0.4207958780492401 0.3596287396468228 0 +-0.436142143265284 -0.2463416916629915 0 +0.3122573415516121 -0.1667305077874104 0 +0.1832065850970164 0.09177111279705727 0 +0.431350132729295 -0.2463969193200146 0 +0.4245259609622692 0.326462497020875 0 +0.3259473550134196 0.4189160364122899 0 +-0.3759437453003407 0.07263414346165893 0 +-0.04999999999875714 0.4351415221593167 0 +0.08902438163851666 0.1709756183613533 0 +-0.2006451802366052 0.3498786945652701 0 +0.2669600775788116 0.3535344791336192 0 +0.3554292542818696 0.2758609329972636 0 +-0.4174064979004342 -0.07317006176782903 0 +0.426794919242815 -0.4267949192429636 0 +-0.4267949192429636 -0.426794919242815 0 +0.09279973326616217 -0.1902773314470203 0 +-0.08684443665693881 -0.183157778505487 0 +-0.4323995522689791 0.2242622428075438 0 +0.07831170102774883 0.03735879558561638 0 +-0.4331911662780547 0.4331911662782387 0 +0.4500000000001133 0.4499999999998865 0 +-0.3553674832650583 -0.1296324926024146 0 +-0.2144365448316186 -0.1744415757602074 0 +-0.4403211998925589 0.05933693233310566 0 +-0.3652998276883476 0.2000898537205348 0 +0.1973037771905921 0.3698060183794389 0 +$EndNodes +$Elements +6 319 1 319 +1 1 1 10 +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 2 +1 2 1 10 +11 2 16 +12 16 17 +13 17 18 +14 18 19 +15 19 20 +16 20 21 +17 21 22 +18 22 23 +19 23 24 +20 24 3 +1 3 1 10 +21 3 25 +22 25 26 +23 26 27 +24 27 28 +25 28 29 +26 29 30 +27 30 31 +28 31 32 +29 32 33 +30 33 4 +1 4 1 10 +31 4 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 1 +1 5 1 5 +41 5 43 +42 43 44 +43 44 45 +44 45 46 +45 46 6 +2 1 2 274 +46 58 107 109 +47 45 44 94 +48 43 44 96 +49 58 109 133 +50 77 113 124 +51 94 44 127 +52 97 56 132 +53 98 97 132 +54 96 44 122 +55 16 17 115 +56 90 55 91 +57 16 115 146 +58 91 51 123 +59 55 90 92 +60 124 113 134 +61 47 109 154 +62 66 5 151 +63 107 58 128 +64 5 66 135 +65 6 25 138 +66 51 91 140 +67 28 55 130 +68 93 50 117 +69 109 110 154 +70 113 77 129 +71 55 92 130 +72 35 98 132 +73 25 26 138 +74 98 35 150 +75 54 61 62 +76 114 50 134 +77 92 90 96 +78 44 43 127 +79 50 93 134 +80 109 47 133 +81 89 145 154 +82 91 55 140 +83 89 95 145 +84 61 59 94 +85 54 21 61 +86 92 96 122 +87 24 6 137 +88 71 104 119 +89 62 61 94 +90 127 43 135 +91 54 63 117 +92 90 87 96 +93 63 93 117 +94 101 112 113 +95 53 105 106 +96 50 114 136 +97 50 19 117 +98 56 97 116 +99 108 106 126 +100 53 9 105 +101 53 106 107 +102 22 59 61 +103 63 62 64 +104 43 96 141 +105 22 23 59 +106 26 27 60 +107 121 49 139 +108 54 62 63 +109 101 111 112 +110 99 100 102 +111 111 57 112 +112 14 57 111 +113 108 126 155 +114 87 83 125 +115 48 99 102 +116 9 10 105 +117 104 102 119 +118 48 11 99 +119 84 83 87 +120 102 100 119 +121 107 106 108 +122 62 94 127 +123 83 80 118 +124 20 21 54 +125 28 29 55 +126 75 73 118 +127 65 64 66 +128 99 52 100 +129 112 57 115 +130 53 107 128 +131 82 80 83 +132 13 14 111 +133 103 102 104 +134 18 19 50 +135 21 22 61 +136 5 125 151 +137 48 102 103 +138 100 52 101 +139 20 54 117 +140 63 64 65 +141 80 75 118 +142 113 112 114 +143 118 73 120 +144 39 40 47 +145 10 11 48 +146 36 37 49 +147 51 31 116 +148 12 52 99 +149 70 68 120 +150 73 70 120 +151 10 48 105 +152 109 108 110 +153 11 12 99 +154 12 13 52 +155 32 56 116 +156 67 66 68 +157 30 31 51 +158 125 120 151 +159 82 83 84 +160 31 32 116 +161 107 108 109 +162 126 76 155 +163 19 20 117 +164 84 87 90 +165 8 9 53 +166 72 70 73 +167 14 15 57 +168 65 66 67 +169 101 113 129 +170 32 33 56 +171 41 42 58 +172 118 120 125 +173 52 13 111 +174 105 103 106 +175 49 121 157 +176 74 73 75 +177 79 75 80 +178 72 73 74 +179 95 38 145 +180 48 103 105 +181 101 52 111 +182 69 68 70 +183 87 125 141 +184 81 80 82 +185 119 100 129 +186 67 68 69 +187 114 112 115 +188 78 86 89 +189 88 97 98 +190 79 80 81 +191 78 75 79 +192 71 70 72 +193 81 85 88 +194 81 82 85 +195 74 75 78 +196 88 85 97 +197 72 74 76 +198 63 65 93 +199 78 79 86 +200 82 84 123 +201 69 70 71 +202 89 86 95 +203 67 69 77 +204 84 91 123 +205 106 103 126 +206 84 90 91 +207 115 57 146 +208 7 8 128 +209 44 45 122 +210 79 81 121 +211 23 24 137 +212 86 79 121 +213 81 88 121 +214 85 82 123 +215 113 114 134 +216 83 118 125 +217 64 62 127 +218 125 5 141 +219 103 104 126 +220 89 110 131 +221 65 67 124 +222 47 40 133 +223 137 46 144 +224 27 28 130 +225 55 29 140 +226 67 77 124 +227 93 65 124 +228 41 58 133 +229 45 94 144 +230 8 53 128 +231 104 71 149 +232 30 51 140 +233 34 35 132 +234 68 66 151 +235 49 37 156 +236 6 46 137 +237 78 89 131 +238 64 127 135 +239 76 74 131 +240 38 39 145 +241 100 101 129 +242 35 36 150 +243 26 60 138 +244 15 2 146 +245 2 16 146 +246 42 1 147 +247 1 7 147 +248 18 50 136 +249 96 87 141 +250 116 97 142 +251 60 27 130 +252 115 17 136 +253 97 85 142 +254 43 5 135 +255 71 119 148 +256 86 121 139 +257 71 72 149 +258 38 95 156 +259 94 59 144 +260 69 71 148 +261 139 49 156 +262 74 78 131 +263 123 51 142 +264 122 45 143 +265 59 23 137 +266 119 129 148 +267 40 41 133 +268 59 137 144 +269 17 18 136 +270 66 64 135 +271 46 6 138 +272 33 4 152 +273 4 34 152 +274 25 6 153 +275 6 24 153 +276 114 115 136 +277 51 116 142 +278 29 30 140 +279 110 89 154 +280 95 86 139 +281 93 124 134 +282 5 43 141 +283 121 88 157 +284 56 33 152 +285 7 128 147 +286 39 47 145 +287 120 68 151 +288 85 123 142 +289 46 138 143 +290 76 126 149 +291 126 104 149 +292 57 15 146 +293 58 42 147 +294 36 49 150 +295 88 98 157 +296 45 46 143 +297 46 45 144 +298 129 77 148 +299 110 108 155 +300 92 122 158 +301 72 76 149 +302 37 38 156 +303 77 69 148 +304 24 3 153 +305 3 25 153 +306 128 58 147 +307 76 131 155 +308 122 143 158 +309 138 60 143 +310 130 92 158 +311 132 56 152 +312 145 47 154 +313 131 110 155 +314 60 130 158 +315 34 132 152 +316 95 139 156 +317 98 150 157 +318 143 60 158 +319 150 49 157 +$EndElements diff --git a/skfem/assembly/__init__.py b/skfem/assembly/__init__.py index beb9da363..975b35379 100644 --- a/skfem/assembly/__init__.py +++ b/skfem/assembly/__init__.py @@ -14,7 +14,7 @@ Named boundaries [# facets]: left [1], bottom [1], right [1], top [1] 2. Create :class:`~skfem.assembly.CellBasis` or - :class:`~skfem.assembly.BoundaryFacetBasis` objects. + :class:`~skfem.assembly.FacetBasis` objects. >>> basis = fem.CellBasis(m, e) diff --git a/skfem/assembly/basis/__init__.py b/skfem/assembly/basis/__init__.py index 4141bf93d..7e6050d3a 100644 --- a/skfem/assembly/basis/__init__.py +++ b/skfem/assembly/basis/__init__.py @@ -1,11 +1,11 @@ from .abstract_basis import AbstractBasis # noqa from .cell_basis import CellBasis # noqa -from .boundary_facet_basis import BoundaryFacetBasis # noqa +from .facet_basis import FacetBasis # noqa from .interior_facet_basis import InteriorFacetBasis # noqa from .mortar_facet_basis import MortarFacetBasis # noqa # aliases Basis = CellBasis InteriorBasis = CellBasis # backwards compatibility -ExteriorFacetBasis = BoundaryFacetBasis # backwards compatibility -FacetBasis = BoundaryFacetBasis +ExteriorFacetBasis = FacetBasis # backwards compatibility +BoundaryFacetBasis = FacetBasis # backwards compatibility diff --git a/skfem/assembly/basis/abstract_basis.py b/skfem/assembly/basis/abstract_basis.py index e63f90998..be0fe4b46 100644 --- a/skfem/assembly/basis/abstract_basis.py +++ b/skfem/assembly/basis/abstract_basis.py @@ -23,7 +23,7 @@ class AbstractBasis: Please see the following implementations: - :class:`~skfem.assembly.CellBasis`, basis functions inside elements - - :class:`~skfem.assembly.ExteriorFacetBasis`, basis functions on boundary + - :class:`~skfem.assembly.FacetBasis`, basis functions on boundary - :class:`~skfem.assembly.InteriorFacetBasis`, basis functions on facets inside the domain @@ -32,12 +32,12 @@ class AbstractBasis: mesh: Mesh elem: Element tind: Optional[ndarray] = None + tind_normals: Optional[ndarray] = None dx: ndarray basis: List[Tuple[DiscreteField, ...]] = [] X: ndarray W: ndarray dofs: Dofs - _side: int = 1 def __init__(self, mesh: Mesh, @@ -240,52 +240,12 @@ def to_indices(f): raise ValueError if elements is not None: - elements = self._normalize_elements(elements) + elements = self.mesh.normalize_elements(elements) return self.dofs.get_element_dofs(elements, skip_dofnames=skip) - facets = self._normalize_facets(facets) + facets = self.mesh.normalize_facets(facets) return self.dofs.get_facet_dofs(facets, skip_dofnames=skip) - def _normalize_facets(self, facets): - if isinstance(facets, ndarray): - return facets - if facets is None: - return self.mesh.boundary_facets() - elif isinstance(facets, (tuple, list, set)): - return np.unique( - np.concatenate( - [self._normalize_facets(f) for f in facets] - ) - ) - elif callable(facets): - return self.mesh.facets_satisfying(facets) - elif isinstance(facets, str): - if ((self.mesh.boundaries is not None - and facets in self.mesh.boundaries)): - return self.mesh.boundaries[facets] - else: - raise ValueError("Boundary '{}' not found.".format(facets)) - raise NotImplementedError - - def _normalize_elements(self, elements): - if isinstance(elements, ndarray): - return elements - if callable(elements): - return self.mesh.elements_satisfying(elements) - elif isinstance(elements, (tuple, list, set)): - return np.unique( - np.concatenate( - [self._normalize_elements(e) for e in elements] - ) - ) - elif isinstance(elements, str): - if ((self.mesh.subdomains is not None - and elements in self.mesh.subdomains)): - return self.mesh.subdomains[elements] - else: - raise ValueError("Subdomain '{}' not found.".format(elements)) - raise NotImplementedError - def __repr__(self): size = sum([sum([y.size if hasattr(y, 'size') else 0 for y in x]) @@ -323,12 +283,17 @@ def interpolate(self, w: ndarray) -> Union[DiscreteField, if w.shape[0] != self.N: raise ValueError("Input array has wrong size.") - refs = self.basis[0] + if isinstance(self.elem, ElementVector): + # ElementVector shouldn't get split here: workaround + refs: List[Tuple[ndarray, 'AbstractBasis']] = [(np.array([]), + self)] + else: + refs = self.split(w) dfs: List[DiscreteField] = [] # loop over solution components for c in range(len(refs)): - ref = refs[c] + ref = refs[c][1].basis[0][0] if ref.is_zero(): dfs.append(ref) continue @@ -384,7 +349,7 @@ def split_indices(self) -> List[ndarray]: e.facet_dofs, e.interior_dofs]) return output - raise ValueError("Basis.elem has only a single component!") + return [np.unique(self.dofs.element_dofs)] def split_bases(self) -> List['AbstractBasis']: """Return Basis objects for the solution components.""" @@ -396,7 +361,7 @@ def split_bases(self) -> List['AbstractBasis']: return [type(self)(self.mesh, self.elem.elem, self.mapping, quadrature=self.quadrature) for _ in range(self.mesh.dim())] - raise ValueError("AbstractBasis.elem has only a single component!") + return [self] @property def quadrature(self): diff --git a/skfem/assembly/basis/cell_basis.py b/skfem/assembly/basis/cell_basis.py index def8766ed..a67b5a4b3 100644 --- a/skfem/assembly/basis/cell_basis.py +++ b/skfem/assembly/basis/cell_basis.py @@ -1,5 +1,5 @@ import logging -from typing import Callable, Optional, Tuple +from typing import Callable, Optional, Tuple, Any import numpy as np from numpy import ndarray @@ -38,7 +38,7 @@ def __init__(self, elem: Element, mapping: Optional[Mapping] = None, intorder: Optional[int] = None, - elements: Optional[ndarray] = None, + elements: Optional[Any] = None, quadrature: Optional[Tuple[ndarray, ndarray]] = None, dofs: Optional[Dofs] = None): """Combine :class:`~skfem.mesh.Mesh` and :class:`~skfem.element.Element` @@ -68,25 +68,27 @@ def __init__(self, logger.info("Initializing {}({}, {})".format(type(self).__name__, type(mesh).__name__, type(elem).__name__)) - super(CellBasis, self).__init__(mesh, - elem, - mapping, - intorder, - quadrature, - mesh.refdom, - dofs) - - self.basis = [self.elem.gbasis(self.mapping, self.X, j, tind=elements) - for j in range(self.Nbfun)] + super(CellBasis, self).__init__( + mesh, + elem, + mapping, + intorder, + quadrature, + mesh.refdom, + dofs, + ) if elements is None: - self.nelems = mesh.nelements self.tind = None + self.nelems = mesh.nelements else: - self.nelems = len(elements) - self.tind = self._normalize_elements(elements) + self.tind = mesh.normalize_elements(elements) + self.nelems = len(self.tind) + + self.basis = [self.elem.gbasis(self.mapping, self.X, j, tind=self.tind) + for j in range(self.Nbfun)] - self.dx = (np.abs(self.mapping.detDF(self.X, tind=elements)) + self.dx = (np.abs(self.mapping.detDF(self.X, tind=self.tind)) * np.tile(self.W, (self.nelems, 1))) logger.info("Initializing finished.") diff --git a/skfem/assembly/basis/boundary_facet_basis.py b/skfem/assembly/basis/facet_basis.py similarity index 83% rename from skfem/assembly/basis/boundary_facet_basis.py rename to skfem/assembly/basis/facet_basis.py index 986cff2ea..ae1592993 100644 --- a/skfem/assembly/basis/boundary_facet_basis.py +++ b/skfem/assembly/basis/facet_basis.py @@ -1,5 +1,5 @@ import logging -from typing import Callable, Optional, Tuple +from typing import Callable, Optional, Tuple, Any import numpy as np from numpy import ndarray @@ -8,6 +8,7 @@ ElementTriP0) from skfem.mapping import Mapping from skfem.mesh import Mesh, MeshHex, MeshLine, MeshQuad, MeshTet, MeshTri +from skfem.generic_utils import OrientedBoundary from .abstract_basis import AbstractBasis from .cell_basis import CellBasis @@ -17,8 +18,8 @@ logger = logging.getLogger(__name__) -class BoundaryFacetBasis(AbstractBasis): - """For fields defined on the boundary of the domain.""" +class FacetBasis(AbstractBasis): + """For integrating over facets of the mesh. Usually over the boundary.""" def __init__(self, mesh: Mesh, @@ -26,9 +27,9 @@ def __init__(self, mapping: Optional[Mapping] = None, intorder: Optional[int] = None, quadrature: Optional[Tuple[ndarray, ndarray]] = None, - facets: Optional[ndarray] = None, + facets: Optional[Any] = None, dofs: Optional[Dofs] = None, - _tind: Optional[ndarray] = None): + side: int = 0): """Precomputed global basis on boundary facets. Parameters @@ -55,39 +56,45 @@ def __init__(self, typestr = ("{}({}, {})".format(type(self).__name__, type(mesh).__name__, type(elem).__name__)) - logger.info("Initializing " + typestr) - super(BoundaryFacetBasis, self).__init__(mesh, - elem, - mapping, - intorder, - quadrature, - mesh.brefdom, - dofs) - - # facets where the basis is evaluated + logger.info("Initializing {}".format(typestr)) + super(FacetBasis, self).__init__( + mesh, + elem, + mapping, + intorder, + quadrature, + mesh.brefdom, + dofs, + ) + + # by default use boundary facets if facets is None: self.find = np.nonzero(self.mesh.f2t[1] == -1)[0] else: - self.find = self._normalize_facets(facets) - - if len(self.find) == 0: - logger.warning("Initializing {} with zero facets.".format(typestr)) + self.find = mesh.normalize_facets(facets) - if _tind is None: - self.tind = self.mesh.f2t[0, self.find] + # fix the orientation + if isinstance(self.find, OrientedBoundary): + self.tind = self.mesh.f2t[(-1) ** side * self.find.ori - side, + self.find] + self.tind_normals = self.mesh.f2t[self.find.ori, self.find] else: - self.tind = _tind + self.tind = self.mesh.f2t[side, self.find] + self.tind_normals = self.mesh.f2t[0, self.find] + + if len(self.find) == 0: + logger.warning("Initializing {} with no facets.".format(typestr)) # boundary refdom to global facet x = self.mapping.G(self.X, find=self.find) # global facet to refdom facet Y = self.mapping.invF(x, tind=self.tind) - # construct normal vectors from side=0 always - Y0 = self.mapping.invF(x, tind=self.mesh.f2t[0, self.find]) + # calculate normals + Y0 = self.mapping.invF(x, tind=self.tind_normals) self.normals = DiscreteField( value=self.mapping.normals(Y0, - self.mesh.f2t[0, self.find], + self.tind_normals, self.find, self.mesh.t2f) ) @@ -120,14 +127,12 @@ def mesh_parameters(self) -> DiscreteField: def _trace_project(self, x: ndarray, elem: Element) -> ndarray: - """Trace mesh basis projection.""" - from skfem.utils import projection - fbasis = BoundaryFacetBasis(self.mesh, - elem, - facets=self.find, - quadrature=(self.X, self.W)) + fbasis = FacetBasis(self.mesh, + elem, + facets=self.find, + quadrature=(self.X, self.W)) I = fbasis.get_dofs(self.find).all() if len(I) == 0: # special case: no facet DOFs if fbasis.dofs.interior_dofs is not None: @@ -205,7 +210,7 @@ def trace(self, self._trace_project(x, target_elem) ) - def with_element(self, elem: Element) -> 'BoundaryFacetBasis': + def with_element(self, elem: Element) -> 'FacetBasis': """Return a similar basis using a different element.""" return type(self)( self.mesh, diff --git a/skfem/assembly/basis/interior_facet_basis.py b/skfem/assembly/basis/interior_facet_basis.py index 307a2101c..d25a49b0a 100644 --- a/skfem/assembly/basis/interior_facet_basis.py +++ b/skfem/assembly/basis/interior_facet_basis.py @@ -1,4 +1,4 @@ -from typing import Optional, Tuple +from typing import Optional, Tuple, Any import numpy as np from numpy import ndarray @@ -6,11 +6,11 @@ from skfem.mapping import Mapping from skfem.mesh import Mesh -from .boundary_facet_basis import BoundaryFacetBasis +from .facet_basis import FacetBasis from ..dofs import Dofs -class InteriorFacetBasis(BoundaryFacetBasis): +class InteriorFacetBasis(FacetBasis): """For evaluating integrals over interior facets. Useful for, e.g., a posteriori error estimators or implementing interior @@ -23,25 +23,23 @@ def __init__(self, mapping: Optional[Mapping] = None, intorder: Optional[int] = None, quadrature: Optional[Tuple[ndarray, ndarray]] = None, - facets: Optional[ndarray] = None, - side: int = 0, - dofs: Optional[Dofs] = None): + facets: Optional[Any] = None, + dofs: Optional[Dofs] = None, + side: int = 0): """Precomputed global basis on interior facets.""" - if side not in (0, 1): - raise Exception("'side' must be 0 or 1.") - if facets is None: facets = np.nonzero(mesh.f2t[1] != -1)[0] - facets = self._normalize_facets(facets) - tind = mesh.f2t[side, facets] - - super(InteriorFacetBasis, self).__init__(mesh, - elem, - mapping=mapping, - intorder=intorder, - quadrature=quadrature, - facets=facets, - _tind=tind, - dofs=dofs) + facets = mesh.normalize_facets(facets) + + super(InteriorFacetBasis, self).__init__( + mesh, + elem, + mapping=mapping, + intorder=intorder, + quadrature=quadrature, + facets=facets, + dofs=dofs, + side=side, + ) diff --git a/skfem/assembly/basis/mortar_facet_basis.py b/skfem/assembly/basis/mortar_facet_basis.py index 95b8b2c51..20c97f49c 100644 --- a/skfem/assembly/basis/mortar_facet_basis.py +++ b/skfem/assembly/basis/mortar_facet_basis.py @@ -1,15 +1,15 @@ -from typing import Optional, Tuple +from typing import Optional, Tuple, Any from numpy import ndarray from skfem.element import Element from skfem.mapping import MappingMortar from skfem.mesh import Mesh -from .boundary_facet_basis import BoundaryFacetBasis +from .facet_basis import FacetBasis from ..dofs import Dofs -class MortarFacetBasis(BoundaryFacetBasis): +class MortarFacetBasis(FacetBasis): def __init__(self, mesh: Mesh, @@ -17,9 +17,9 @@ def __init__(self, mapping: MappingMortar, intorder: Optional[int] = None, quadrature: Optional[Tuple[ndarray, ndarray]] = None, - facets: Optional[ndarray] = None, - side: int = 0, - dofs: Optional[Dofs] = None): + facets: Optional[Any] = None, + dofs: Optional[Dofs] = None, + side: int = 0): """Precomputed global basis on the mortar mesh.""" if side not in (0, 1): @@ -31,12 +31,14 @@ def __init__(self, mapping.side = side facets = mapping.helper_to_orig[side] - facets = self._normalize_facets(facets) - - super(MortarFacetBasis, self).__init__(mesh, - elem, - mapping=mapping, - intorder=intorder, - quadrature=quadrature, - facets=facets, - dofs=dofs) + facets = mesh.normalize_facets(facets) + + super(MortarFacetBasis, self).__init__( + mesh, + elem, + mapping=mapping, + intorder=intorder, + quadrature=quadrature, + facets=facets, + dofs=dofs, + ) diff --git a/skfem/element/__init__.py b/skfem/element/__init__.py index da83cb74d..bb0f57ea8 100644 --- a/skfem/element/__init__.py +++ b/skfem/element/__init__.py @@ -3,7 +3,7 @@ In order to use an element, you simply initialize the respective object and pass it to the constructor of :class:`~skfem.assembly.CellBasis` or -:class:`~skfem.assembly.BoundaryFacetBasis`. See below for a list of supported +:class:`~skfem.assembly.FacetBasis`. See below for a list of supported elements. """ diff --git a/skfem/generic_utils.py b/skfem/generic_utils.py index dee3c9d8f..f74c825c4 100644 --- a/skfem/generic_utils.py +++ b/skfem/generic_utils.py @@ -1,3 +1,5 @@ +import numpy as np + from numpy import ndarray @@ -6,3 +8,18 @@ def hash_args(*args): return tuple(hash(arg.tobytes()) if isinstance(arg, ndarray) else hash(arg) for arg in args) + + +class OrientedBoundary(ndarray): + """An array of facet indices with orientation.""" + + def __new__(cls, indices, ori): + obj = np.asarray(indices).view(cls) + obj.ori = np.array(ori, dtype=int) + assert len(obj) == len(obj.ori) + return obj + + def __array_finalize__(self, obj): + if obj is None: + return + self.ori = getattr(obj, 'ori', None) diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py index d2a913b64..0fa894dfd 100644 --- a/skfem/io/meshio.py +++ b/skfem/io/meshio.py @@ -2,20 +2,23 @@ import meshio import numpy as np -import skfem + +from skfem.mesh import (MeshTet1, MeshTet2, MeshHex1, MeshHex2, MeshWedge1, + MeshTri1, MeshTri2, MeshQuad1, MeshQuad2, MeshLine1) +from skfem.generic_utils import OrientedBoundary MESH_TYPE_MAPPING = { - 'tetra': skfem.MeshTet1, - 'tetra10': skfem.MeshTet2, - 'hexahedron': skfem.MeshHex1, - 'hexahedron27': skfem.MeshHex2, - 'wedge': skfem.MeshWedge1, - 'triangle': skfem.MeshTri1, - 'triangle6': skfem.MeshTri2, - 'quad': skfem.MeshQuad1, - 'quad9': skfem.MeshQuad2, - 'line': skfem.MeshLine1, + 'tetra': MeshTet1, + 'tetra10': MeshTet2, + 'hexahedron': MeshHex1, + 'hexahedron27': MeshHex2, + 'wedge': MeshWedge1, + 'triangle': MeshTri1, + 'triangle6': MeshTri2, + 'quad': MeshQuad1, + 'quad9': MeshQuad2, + 'line': MeshLine1, } BOUNDARY_TYPE_MAPPING = { @@ -114,18 +117,47 @@ def from_meshio(m, # parse boundaries from cell_sets if m.cell_sets and bnd_type in m.cells_dict: - facets = { - k: [tuple(f) for f in np.sort(m.cells_dict[bnd_type][v[bnd_type]])] + oriented_facets = { + k: [tuple(f) for f in m.cells_dict[bnd_type][v[bnd_type]]] for k, v in m.cell_sets_dict.items() if bnd_type in v and k.split(":")[0] != "gmsh" } - boundaries = {k: np.array([i for i, f in - enumerate(map(tuple, mtmp.facets.T)) - if f in v]) - for k, v in facets.items()} + sorted_facets = {k: [tuple(np.sort(f)) for f in v] + for k, v in oriented_facets.items()} + + for k, v in oriented_facets.items(): + indices = [] + oris = [] + for i, f in enumerate(map(tuple, mtmp.facets.T)): + if f in sorted_facets[k]: + indices.append(i) + ix = sorted_facets[k].index(f) + facet = v[ix] + t1, t2 = mtmp.f2t[:, i] + if t2 == -1: + # orientation on boundary is 0 + oris.append(0) + continue + if len(f) == 2: + # rotate tangent to find normal + tangent = mtmp.p[:, facet[1]] - mtmp.p[:, facet[0]] + normal = np.array([-tangent[1], tangent[0]]) + elif len(f) == 3: + # cross product to find normal + tangent1 = mtmp.p[:, facet[1]] - mtmp.p[:, facet[0]] + tangent2 = mtmp.p[:, facet[2]] - mtmp.p[:, facet[0]] + normal = -np.cross(tangent1, tangent2) + else: + raise NotImplementedError + # find another vector using node outside of boundary + third = np.setdiff1d(mtmp.t[:, t1], + np.array(f))[0] + outplane = mtmp.p[:, f[0]] - mtmp.p[:, third] + oris.append(1 if np.dot(normal, outplane) > 0 else 0) + boundaries[k] = OrientedBoundary(indices, oris) # MSH 2.2 tag parsing - if m.cell_data and m.field_data: + if len(boundaries) == 0 and m.cell_data and m.field_data: try: elements_tag = m.cell_data_dict['gmsh:physical'][meshio_type] subdomains = {} @@ -196,9 +228,9 @@ def to_meshio(mesh, encode_point_data=False): t = mesh.dofs.element_dofs.copy() - if isinstance(mesh, skfem.MeshHex2): + if isinstance(mesh, MeshHex2): t = t[HEX_MAPPING] - elif isinstance(mesh, skfem.MeshHex): + elif isinstance(mesh, MeshHex1): t = t[HEX_MAPPING[:8]] mtype = TYPE_MESH_MAPPING[type(mesh)] @@ -232,6 +264,9 @@ def to_file(mesh, encode_point_data=False, **kwargs): + if 'file_format' not in kwargs and filename.split('.')[-1] == 'msh': + kwargs.update({'file_format': 'gmsh'}) + meshio.write(filename, to_meshio(mesh, point_data, diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py index 0a9e7d44d..6fd679018 100644 --- a/skfem/mesh/mesh.py +++ b/skfem/mesh/mesh.py @@ -8,6 +8,7 @@ from numpy import ndarray from ..element import BOUNDARY_ELEMENT_MAP, Element +from ..generic_utils import OrientedBoundary logger = logging.getLogger(__name__) @@ -146,16 +147,19 @@ def boundary_edges(self) -> ndarray: return edge_candidates[ix] def with_boundaries(self, - boundaries: Dict[str, Callable[[ndarray], ndarray]], + boundaries: Dict[str, Union[Callable[[ndarray], + ndarray], + ndarray]], boundaries_only: bool = True): """Return a copy of the mesh with named boundaries. Parameters ---------- boundaries - A dictionary of lambda functions with the names of the boundaries - as keys. The midpoint of the facet should return ``True`` for the - corresponding lambda function if the facet belongs to the boundary. + A dictionary of index arrays or lambda functions with the names of + the boundaries as keys. The midpoint of the facet should return + ``True`` for the corresponding lambda function if the facet belongs + to the boundary. boundaries_only If ``True``, consider only facets on the boundary of the domain. @@ -164,8 +168,9 @@ def with_boundaries(self, self, _boundaries={ **({} if self._boundaries is None else self._boundaries), - **{name: self.facets_satisfying(test, boundaries_only) - for name, test in boundaries.items()} + **{name: self.facets_satisfying(test_or_set, boundaries_only) + if callable(test_or_set) else test_or_set + for name, test_or_set in boundaries.items()} }, ) @@ -284,7 +289,8 @@ def nodes_satisfying(self, def facets_satisfying(self, test: Callable[[ndarray], ndarray], - boundaries_only: bool = False) -> ndarray: + boundaries_only: bool = False, + normal: Optional[ndarray] = None) -> ndarray: """Return facets whose midpoints satisfy some condition. Parameters @@ -294,14 +300,46 @@ def facets_satisfying(self, to be included in the return set. boundaries_only If ``True``, include only boundary facets. + normal + If given, used to orient the set of facets. """ midp = self.p[:, self.facets].mean(axis=1) facets = np.nonzero(test(midp))[0] if boundaries_only: facets = np.intersect1d(facets, self.boundary_facets()) + if normal is not None: + tind = self.f2t[0, facets] + mapping = self._mapping() + normals = mapping.normals(np.zeros((self.dim(), 1)), + tind, + facets, + self.t2f).T[0].T + ori = 1 * (np.dot(normal, normals) < 0) + return OrientedBoundary(facets, ori) return facets + def facets_around(self, elements, flip=False) -> OrientedBoundary: + """Return the oriented set of facets around a set of elements. + + Parameters + ---------- + elements + An array of element indices or, alternatively, something that can + be cast into one via ``Mesh.normalize_elements``. + flip + If ``True``, use traces outside the subdomain and inward normals. + + """ + elements = self.normalize_elements(elements) + facets, counts = np.unique(self.t2f[:, elements], return_counts=True) + facets = facets[counts == 1] + if flip: + ori = np.nonzero(~np.isin(self.f2t[:, facets], elements).T)[1] + else: + ori = np.nonzero(np.isin(self.f2t[:, facets], elements).T)[1] + return OrientedBoundary(facets, ori) + def elements_satisfying(self, test: Callable[[ndarray], ndarray]) -> ndarray: """Return elements whose midpoints satisfy some condition. @@ -921,3 +959,80 @@ def plot(self, x, visuals='matplotlib', **kwargs): """Convenience wrapper for skfem.visuals.""" mod = importlib.import_module('skfem.visuals.{}'.format(visuals)) return mod.plot(self, x, **kwargs) + + def normalize_facets(self, facets) -> ndarray: + """Generate an array of facet indices. + + Parameters + ---------- + elements + Criteria of which facets to include. Function has different + behavior based on the type of this parameter. + + """ + if isinstance(facets, int): + # Make normalize_facets([1,2,3]) have the same behavior as + # normalize_facets(np.array([1,2,3])) + return np.array([facets]) + if isinstance(facets, ndarray): + # Assume the facets have already been normalized + return facets + if facets is None: + # Default behavior. + return self.boundary_facets() + elif isinstance(facets, (tuple, list, set)): + # Recurse over the list, building an array of all matching facets + return np.unique( + np.concatenate( + [self.normalize_facets(f) for f in facets] + ) + ) + elif callable(facets): + # The callable should accept an array of facet centers and return + # an boolean array with True for facets that should be included. + return self.facets_satisfying(facets) + elif isinstance(facets, str): + # Assume string is the label of a boundary in the mesh. + if ((self.boundaries is not None + and facets in self.boundaries)): + return self.boundaries[facets] + else: + raise ValueError("Boundary '{}' not found.".format(facets)) + raise NotImplementedError + + def normalize_elements(self, elements) -> ndarray: + """Generate an array of element indices. + + Parameters + ---------- + elements + Criteria of which elements to include. Function has different + behavior based on the type of this parameter. + + """ + if isinstance(elements, int): + # Make normalize_elements([1,2,3]) have the same behavior as + # normalize_elements(np.array([1,2,3])) + return np.array([elements]) + if isinstance(elements, ndarray): + # Assume the elements have already been normalized + return elements + if callable(elements): + # The callable should accept an array of element centers and return + # an boolean array with True for elements that should be included. + return self.elements_satisfying(elements) + elif isinstance(elements, (tuple, list, set)): + # Recurse over the list, building an array of all matching elements + return np.unique( + np.concatenate( + [self.normalize_elements(e) for e in elements] + ) + ) + elif isinstance(elements, str): + # Assume string is the label of a subdomain in the mesh. + if ((self.subdomains is not None + and elements in self.subdomains)): + return self.subdomains[elements] + else: + raise ValueError("Subdomain '{}' not found.".format(elements)) + raise NotImplementedError diff --git a/skfem/visuals/matplotlib.py b/skfem/visuals/matplotlib.py index 9d95e3b79..2c018b34f 100644 --- a/skfem/visuals/matplotlib.py +++ b/skfem/visuals/matplotlib.py @@ -7,6 +7,7 @@ import matplotlib.pyplot as plt from matplotlib.axes import Axes +from matplotlib.collections import PolyCollection from ..assembly import CellBasis from ..mesh import Mesh2D, MeshLine1, MeshQuad1, MeshTri1, Mesh3D @@ -129,7 +130,13 @@ def draw_mesh2d(m: Mesh2D, **kwargs) -> Axes: ys, color=colors[i % len(colors)], linewidth=(kwargs['linewidth'] - if 'linewidth' in kwargs else 1.)) + if 'linewidth' in kwargs else 2.)) + if hasattr(m.boundaries[k], 'ori'): + tris = m.f2t[m.boundaries[k].ori, m.boundaries[k]] + color = colors[i % len(colors)][:3] + (.1,) + collec = PolyCollection(m.p[:, m.t[:, tris]].T, + color=color) + ax.add_collection(collec) if "node_numbering" in kwargs: for itr in range(m.p.shape[1]): diff --git a/tests/test_basis.py b/tests/test_basis.py index 8bef172ca..14f47ce1d 100644 --- a/tests/test_basis.py +++ b/tests/test_basis.py @@ -1,5 +1,6 @@ import pickle from unittest import TestCase +from pathlib import Path import pytest import numpy as np @@ -17,10 +18,17 @@ ElementLineP0, ElementQuad1, ElementQuad0, ElementTetP1, ElementTetP0, ElementHex1, ElementHex0, ElementLineP1, ElementLineMini, - ElementWedge1, ElementTriRT0, ElementQuadRT0) + ElementWedge1, ElementTriRT0, ElementQuadRT0, + ElementTriP1) +from skfem.helpers import dot, grad +from skfem.utils import enforce +from skfem.io.meshio import to_meshio, from_meshio from skfem.models.poisson import laplace +MESH_PATH = Path(__file__).parents[1] / 'docs' / 'examples' / 'meshes' + + class TestCompositeSplitting(TestCase): def runTest(self): @@ -60,6 +68,7 @@ def bilinf(u, p, v, q, w): U, P = basis.interpolate(x) self.assertTrue(isinstance(U.value, np.ndarray)) self.assertTrue(isinstance(P.value, np.ndarray)) + self.assertTrue(P.shape[0] == m.nelements) self.assertTrue((basis.doflocs[:, D['up'].all()][1] == 1.).all()) @@ -442,3 +451,168 @@ def int_y(w): return w.y ** 2 assert_almost_equal(int_y.assemble(basis, y=y), 1.) + + +def test_subdomain_facet_assembly(): + + def subdomain(x): + return np.logical_and( + np.logical_and(x[0]>.25, x[0]<.75), + np.logical_and(x[1]>.25, x[1]<.75), + ) + + m, e = MeshTri().refined(4), ElementTriP2() + cbasis = CellBasis(m, e) + cbasis_p0 = cbasis.with_element(ElementTriP0()) + + sfbasis = FacetBasis(m, e, facets=m.facets_around(subdomain, flip=True)) + sfbasis_p0 = sfbasis.with_element(ElementTriP0()) + sigma = cbasis_p0.zeros() + 1 + + @BilinearForm + def laplace(u, v, w): + return dot(w.sigma * grad(u), grad(v)) + + A = laplace.assemble(cbasis, sigma=cbasis_p0.interpolate(sigma)) + u0 = cbasis.zeros() + u0[cbasis.get_dofs(elements=subdomain)] = 1 + u0_dofs = cbasis.get_dofs() + cbasis.get_dofs(elements=subdomain) + A, b = enforce(A, D=u0_dofs, x=u0) + u = solve(A, b) + + @Functional + def measure_current(w): + return dot(w.n, w.sigma * grad(w.u)) + + meas = measure_current.assemble(sfbasis, + sigma=sfbasis_p0.interpolate(sigma), + u=sfbasis.interpolate(u)) + + assert_almost_equal(meas, 9.751915526759191) + + +def test_subdomain_facet_assembly_2(): + + m = MeshTri().refined(4).with_subdomains({'all': lambda x: x[0] * 0 + 1}) + e = ElementTriP1() + + @Functional + def boundary_integral(w): + return dot(w.n, grad(w.u)) + + sfbasis = FacetBasis(m, e, facets=m.facets_around('all')) + fbasis = FacetBasis(m, e) + + assert_almost_equal( + boundary_integral.assemble(sfbasis, u=m.p[0] * m.p[1]), + boundary_integral.assemble(fbasis, u=m.p[0] * m.p[1]), + ) + + +@pytest.mark.parametrize( + "m, e, facets, fun", + [ + ( + MeshTri.load(MESH_PATH / 'interface.msh'), + ElementTriP1(), + 'interfacee', + lambda m: m.p[1], + ), + ( + MeshTet.load(MESH_PATH / 'cuubat.msh'), + ElementTetP1(), + 'interface', + lambda m: m.p[0], + ) + ] +) +def test_oriented_interface_integral(m, e, facets, fun): + + fb = FacetBasis(m, e, facets=facets) + assert_almost_equal( + Functional(lambda w: dot(w.fun.grad, w.n)).assemble(fb, fun=fun(m)), + 1.0, + ) + + +@pytest.mark.parametrize( + "m, e, facets, normal", + [ + ( + MeshTri().refined(2) + MeshTri().translated((1.0, 0.0)).refined(2), + ElementTriP1(), + lambda x: x[0] == 1.0, + np.array([1, 0]), + ), + ( + MeshHex().refined(2) + MeshHex().translated((1., 0., 0.)).refined(2), + ElementHex1(), + lambda x: x[0] == 1.0, + np.array([1, 0, 0]), + ), + ] +) +def test_oriented_interface_integral2(m, e, facets, normal): + + fb = FacetBasis(m, e, facets=m.facets_satisfying(facets, normal=normal)) + assert_almost_equal( + Functional(lambda w: dot(w.fun.grad, w.n)).assemble(fb, fun=m.p[0]), + 1.0, + ) + + +@pytest.mark.parametrize( + "m, e", + [ + ( + MeshTri().refined(6).with_subdomains({ + 'mid': lambda x: (x[0] - .5) ** 2 + (x[1] - .5) ** 2 < .5 ** 2, + }), + ElementTriP1(), + ), + ( + MeshTet().refined(4).with_subdomains({ + 'mid': lambda x: ((x[0] - .5) ** 2 + + (x[1] - .5) ** 2 + + (x[2] - .5) ** 2) < .5 ** 2, + }), + ElementTetP1(), + ), + ( + MeshHex().refined(4).with_subdomains({ + 'mid': lambda x: ((x[0] - .5) ** 2 + + (x[1] - .5) ** 2 + + (x[2] - .5) ** 2) < .5 ** 2, + }), + ElementHex1(), + ), + ] +) +def test_oriented_gauss_integral(m, e): + + facets = m.facets_around('mid') + fb = FacetBasis(m, e, facets=facets) + cb = CellBasis(m, e, elements='mid') + assert_almost_equal( + Functional(lambda w: w.x[0] * w.n[0]).assemble(fb), + Functional(lambda w: 1. + 0. * w.x[0]).assemble(cb), + decimal=5, + ) + + +## TODO preserve orientation in save/load cycle +# def test_oriented_saveload(): + +# m = MeshTri().refined(4) +# m = m.with_boundaries({ +# 'mid': m.facets_around([5]), +# }) +# assert len(m.boundaries['mid'].ori) == 3 + +# # cycle to meshio and check that orientation is preserved +# M = from_meshio(to_meshio(m)) + +# assert_almost_equal( +# m.boundaries['mid'].ori, +# M.boundaries['mid'].ori, +# ) diff --git a/tests/test_mesh.py b/tests/test_mesh.py index be2333451..41951b197 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -16,6 +16,9 @@ from skfem.io.json import to_dict, from_dict +MESH_PATH = Path(__file__).parents[1] / 'docs' / 'examples' / 'meshes' + + class MeshTests(TestCase): """Test some of the methods in mesh classes that are not tested elsewhere.""" @@ -61,9 +64,6 @@ def _runTest(self): # disabled np.array([[0.0, 1.0, 2.0], [1.0, 2.0, 3.0]]).T) -MESH_PATH = Path(__file__).parents[1] / 'docs' / 'examples' / 'meshes' - - class Loading(TestCase): """Check that Mesh.load works properly.""" @@ -528,9 +528,9 @@ def test_saveload_cycle_vtk(m): from tempfile import NamedTemporaryFile m = m.refined(2) - with NamedTemporaryFile() as f: - m.save(f.name + ".vtk") - m2 = Mesh.load(f.name + ".vtk") + with NamedTemporaryFile(suffix='.vtk') as f: + m.save(f.name) + m2 = Mesh.load(f.name) assert_array_equal(m.p, m2.p) assert_array_equal(m.t, m2.t) @@ -564,10 +564,10 @@ def test_saveload_cycle_tags(fmt, kwargs, m): .with_boundaries({'test': lambda x: (x[0] == 0) * (x[1] < 0.6), 'set': lambda x: (x[0] == 0) * (x[1] > 0.3)})) from tempfile import NamedTemporaryFile - with NamedTemporaryFile() as f: - m.save(f.name + fmt, point_data={'foo': m.p[0]}, **kwargs) + with NamedTemporaryFile(suffix=fmt) as f: + m.save(f.name, point_data={'foo': m.p[0]}, **kwargs) out = ['point_data', 'cells_dict'] - m2 = Mesh.load(f.name + fmt, out=out) + m2 = Mesh.load(f.name, out=out) assert_array_equal(m.p, m2.p)