dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #18016
Merge request for some sandbox stuff
-
To:
DOLFIN Mailing List <dolfin@xxxxxxxxxxxxxxxxxxx>
-
From:
Andre Massing <massing@xxxxxxxxx>
-
Date:
Tue, 6 Apr 2010 11:25:07 +0200
-
Organization:
Simula Research Laboratory
-
User-agent:
KMail/1.13.1 (Linux/2.6.31.4-centrino-duo; KDE/4.4.1; x86_64; ; )
Hi!
While I am working on a implementation of Nitsche methods for overlapping
meshes in DOLFIN I produce some sand, which might be nice to be added to the
sandbox :) It contains among other things some prototypes for building a
complete intersection map for two overlapping meshes (which cells/facets are
covered by which cells/facest) and for barycenter quadrature formulas on
arbitrary polyhedrons and polygons (unfortunately just for the 3D case, up to
now).
Kind regards,
Andre
# Bazaar merge directive format 2 (Bazaar 0.90)
# revision_id: massing@xxxxxxxxx-20100406091417-ogh30d38n95eh5cq
# target_branch: bzr+ssh://bazaar.launchpad.net/~dolfin-\
# core/dolfin/sandbox/
# testament_sha1: 74c2e1fd30a903bacbd8c50606cb683173e73813
# timestamp: 2010-04-06 11:15:17 +0200
# base_revision_id: logg@xxxxxxxxx-20100226124738-bk9jlyb9voycqbwa
#
# Begin patch
=== renamed directory 'cgal' => 'nitsche_method'
=== added directory 'nitsche_method/data'
=== added directory 'nitsche_method/data/meshes'
=== added file 'nitsche_method/data/meshes/rotator.xml'
--- nitsche_method/data/meshes/rotator.xml 1970-01-01 00:00:00 +0000
+++ nitsche_method/data/meshes/rotator.xml 2010-03-16 14:39:09 +0000
@@ -0,0 +1,402 @@
+<?xml version="1.0" encoding="UTF-8"?>
+
+<dolfin xmlns:dolfin="http://www.fenics.org/dolfin/";>
+ <mesh celltype="triangle" dim="2">
+ <vertices size="117">
+ <vertex index="0" x="-0.1" y="-1"/>
+ <vertex index="1" x="0" y="-1"/>
+ <vertex index="2" x="0" y="-0.9"/>
+ <vertex index="3" x="-0.1" y="-0.9"/>
+ <vertex index="4" x="0.1" y="-1"/>
+ <vertex index="5" x="0.1" y="-0.9"/>
+ <vertex index="6" x="0" y="-0.8"/>
+ <vertex index="7" x="-0.1" y="-0.8"/>
+ <vertex index="8" x="0.1" y="-0.8"/>
+ <vertex index="9" x="0" y="-0.7"/>
+ <vertex index="10" x="-0.1" y="-0.7"/>
+ <vertex index="11" x="0.1" y="-0.7"/>
+ <vertex index="12" x="0" y="-0.6"/>
+ <vertex index="13" x="-0.1" y="-0.6"/>
+ <vertex index="14" x="0.1" y="-0.6"/>
+ <vertex index="15" x="0" y="-0.5"/>
+ <vertex index="16" x="-0.1" y="-0.5"/>
+ <vertex index="17" x="0.1" y="-0.5"/>
+ <vertex index="18" x="0" y="-0.4"/>
+ <vertex index="19" x="-0.1" y="-0.4"/>
+ <vertex index="20" x="0.1" y="-0.4"/>
+ <vertex index="21" x="0" y="-0.3"/>
+ <vertex index="22" x="-0.1" y="-0.3"/>
+ <vertex index="23" x="0.1" y="-0.3"/>
+ <vertex index="24" x="0" y="-0.2"/>
+ <vertex index="25" x="-0.1" y="-0.2"/>
+ <vertex index="26" x="0.1" y="-0.2"/>
+ <vertex index="27" x="0" y="-0.1"/>
+ <vertex index="28" x="-0.1" y="-0.1"/>
+ <vertex index="29" x="0.1" y="-0.1"/>
+ <vertex index="30" x="-1" y="-0.1"/>
+ <vertex index="31" x="-0.9" y="-0.1"/>
+ <vertex index="32" x="-0.9" y="0"/>
+ <vertex index="33" x="-1" y="0"/>
+ <vertex index="34" x="-0.8" y="-0.1"/>
+ <vertex index="35" x="-0.8" y="0"/>
+ <vertex index="36" x="-0.7" y="-0.1"/>
+ <vertex index="37" x="-0.7" y="0"/>
+ <vertex index="38" x="-0.6" y="-0.1"/>
+ <vertex index="39" x="-0.6" y="0"/>
+ <vertex index="40" x="-0.5" y="-0.1"/>
+ <vertex index="41" x="-0.5" y="0"/>
+ <vertex index="42" x="-0.4" y="-0.1"/>
+ <vertex index="43" x="-0.4" y="0"/>
+ <vertex index="44" x="-0.3" y="-0.1"/>
+ <vertex index="45" x="-0.3" y="0"/>
+ <vertex index="46" x="-0.2" y="-0.1"/>
+ <vertex index="47" x="-0.2" y="0"/>
+ <vertex index="48" x="-0.1" y="0"/>
+ <vertex index="49" x="0" y="0"/>
+ <vertex index="50" x="0.1" y="0"/>
+ <vertex index="51" x="0.2" y="-0.1"/>
+ <vertex index="52" x="0.2" y="0"/>
+ <vertex index="53" x="0.3" y="-0.1"/>
+ <vertex index="54" x="0.3" y="0"/>
+ <vertex index="55" x="0.4" y="-0.1"/>
+ <vertex index="56" x="0.4" y="0"/>
+ <vertex index="57" x="0.5" y="-0.1"/>
+ <vertex index="58" x="0.5" y="0"/>
+ <vertex index="59" x="0.6" y="-0.1"/>
+ <vertex index="60" x="0.6" y="0"/>
+ <vertex index="61" x="0.7" y="-0.1"/>
+ <vertex index="62" x="0.7" y="0"/>
+ <vertex index="63" x="0.8" y="-0.1"/>
+ <vertex index="64" x="0.8" y="0"/>
+ <vertex index="65" x="0.9" y="-0.1"/>
+ <vertex index="66" x="0.9" y="0"/>
+ <vertex index="67" x="1" y="-0.1"/>
+ <vertex index="68" x="1" y="0"/>
+ <vertex index="69" x="-0.9" y="0.1"/>
+ <vertex index="70" x="-1" y="0.1"/>
+ <vertex index="71" x="-0.8" y="0.1"/>
+ <vertex index="72" x="-0.7" y="0.1"/>
+ <vertex index="73" x="-0.6" y="0.1"/>
+ <vertex index="74" x="-0.5" y="0.1"/>
+ <vertex index="75" x="-0.4" y="0.1"/>
+ <vertex index="76" x="-0.3" y="0.1"/>
+ <vertex index="77" x="-0.2" y="0.1"/>
+ <vertex index="78" x="-0.1" y="0.1"/>
+ <vertex index="79" x="0" y="0.1"/>
+ <vertex index="80" x="0.1" y="0.1"/>
+ <vertex index="81" x="0.2" y="0.1"/>
+ <vertex index="82" x="0.3" y="0.1"/>
+ <vertex index="83" x="0.4" y="0.1"/>
+ <vertex index="84" x="0.5" y="0.1"/>
+ <vertex index="85" x="0.6" y="0.1"/>
+ <vertex index="86" x="0.7" y="0.1"/>
+ <vertex index="87" x="0.8" y="0.1"/>
+ <vertex index="88" x="0.9" y="0.1"/>
+ <vertex index="89" x="1" y="0.1"/>
+ <vertex index="90" x="0" y="0.2"/>
+ <vertex index="91" x="-0.1" y="0.2"/>
+ <vertex index="92" x="0.1" y="0.2"/>
+ <vertex index="93" x="0" y="0.3"/>
+ <vertex index="94" x="-0.1" y="0.3"/>
+ <vertex index="95" x="0.1" y="0.3"/>
+ <vertex index="96" x="0" y="0.4"/>
+ <vertex index="97" x="-0.1" y="0.4"/>
+ <vertex index="98" x="0.1" y="0.4"/>
+ <vertex index="99" x="0" y="0.5"/>
+ <vertex index="100" x="-0.1" y="0.5"/>
+ <vertex index="101" x="0.1" y="0.5"/>
+ <vertex index="102" x="0" y="0.6"/>
+ <vertex index="103" x="-0.1" y="0.6"/>
+ <vertex index="104" x="0.1" y="0.6"/>
+ <vertex index="105" x="0" y="0.7"/>
+ <vertex index="106" x="-0.1" y="0.7"/>
+ <vertex index="107" x="0.1" y="0.7"/>
+ <vertex index="108" x="0" y="0.8"/>
+ <vertex index="109" x="-0.1" y="0.8"/>
+ <vertex index="110" x="0.1" y="0.8"/>
+ <vertex index="111" x="0" y="0.9"/>
+ <vertex index="112" x="-0.1" y="0.9"/>
+ <vertex index="113" x="0.1" y="0.9"/>
+ <vertex index="114" x="0" y="1"/>
+ <vertex index="115" x="-0.1" y="1"/>
+ <vertex index="116" x="0.1" y="1"/>
+ </vertices>
+ <cells size="152">
+ <triangle index="0" v0="0" v1="1" v2="2"/>
+ <triangle index="1" v0="0" v1="2" v2="3"/>
+ <triangle index="2" v0="1" v1="4" v2="5"/>
+ <triangle index="3" v0="1" v1="2" v2="5"/>
+ <triangle index="4" v0="2" v1="3" v2="6"/>
+ <triangle index="5" v0="3" v1="6" v2="7"/>
+ <triangle index="6" v0="2" v1="5" v2="8"/>
+ <triangle index="7" v0="2" v1="6" v2="8"/>
+ <triangle index="8" v0="6" v1="7" v2="9"/>
+ <triangle index="9" v0="7" v1="9" v2="10"/>
+ <triangle index="10" v0="6" v1="8" v2="11"/>
+ <triangle index="11" v0="6" v1="9" v2="11"/>
+ <triangle index="12" v0="9" v1="10" v2="12"/>
+ <triangle index="13" v0="10" v1="12" v2="13"/>
+ <triangle index="14" v0="9" v1="11" v2="14"/>
+ <triangle index="15" v0="9" v1="12" v2="14"/>
+ <triangle index="16" v0="12" v1="13" v2="15"/>
+ <triangle index="17" v0="13" v1="15" v2="16"/>
+ <triangle index="18" v0="12" v1="14" v2="17"/>
+ <triangle index="19" v0="12" v1="15" v2="17"/>
+ <triangle index="20" v0="15" v1="16" v2="18"/>
+ <triangle index="21" v0="16" v1="18" v2="19"/>
+ <triangle index="22" v0="15" v1="17" v2="20"/>
+ <triangle index="23" v0="15" v1="18" v2="20"/>
+ <triangle index="24" v0="18" v1="19" v2="21"/>
+ <triangle index="25" v0="19" v1="21" v2="22"/>
+ <triangle index="26" v0="18" v1="20" v2="23"/>
+ <triangle index="27" v0="18" v1="21" v2="23"/>
+ <triangle index="28" v0="21" v1="22" v2="24"/>
+ <triangle index="29" v0="22" v1="24" v2="25"/>
+ <triangle index="30" v0="21" v1="23" v2="26"/>
+ <triangle index="31" v0="21" v1="24" v2="26"/>
+ <triangle index="32" v0="24" v1="25" v2="27"/>
+ <triangle index="33" v0="25" v1="27" v2="28"/>
+ <triangle index="34" v0="24" v1="26" v2="29"/>
+ <triangle index="35" v0="24" v1="27" v2="29"/>
+ <triangle index="36" v0="30" v1="31" v2="32"/>
+ <triangle index="37" v0="30" v1="32" v2="33"/>
+ <triangle index="38" v0="31" v1="34" v2="35"/>
+ <triangle index="39" v0="31" v1="32" v2="35"/>
+ <triangle index="40" v0="34" v1="36" v2="37"/>
+ <triangle index="41" v0="34" v1="35" v2="37"/>
+ <triangle index="42" v0="36" v1="38" v2="39"/>
+ <triangle index="43" v0="36" v1="37" v2="39"/>
+ <triangle index="44" v0="38" v1="40" v2="41"/>
+ <triangle index="45" v0="38" v1="39" v2="41"/>
+ <triangle index="46" v0="40" v1="42" v2="43"/>
+ <triangle index="47" v0="40" v1="41" v2="43"/>
+ <triangle index="48" v0="42" v1="44" v2="45"/>
+ <triangle index="49" v0="42" v1="43" v2="45"/>
+ <triangle index="50" v0="44" v1="46" v2="47"/>
+ <triangle index="51" v0="44" v1="45" v2="47"/>
+ <triangle index="52" v0="28" v1="46" v2="48"/>
+ <triangle index="53" v0="46" v1="47" v2="48"/>
+ <triangle index="54" v0="27" v1="28" v2="49"/>
+ <triangle index="55" v0="28" v1="48" v2="49"/>
+ <triangle index="56" v0="27" v1="29" v2="50"/>
+ <triangle index="57" v0="27" v1="49" v2="50"/>
+ <triangle index="58" v0="29" v1="51" v2="52"/>
+ <triangle index="59" v0="29" v1="50" v2="52"/>
+ <triangle index="60" v0="51" v1="53" v2="54"/>
+ <triangle index="61" v0="51" v1="52" v2="54"/>
+ <triangle index="62" v0="53" v1="55" v2="56"/>
+ <triangle index="63" v0="53" v1="54" v2="56"/>
+ <triangle index="64" v0="55" v1="57" v2="58"/>
+ <triangle index="65" v0="55" v1="56" v2="58"/>
+ <triangle index="66" v0="57" v1="59" v2="60"/>
+ <triangle index="67" v0="57" v1="58" v2="60"/>
+ <triangle index="68" v0="59" v1="61" v2="62"/>
+ <triangle index="69" v0="59" v1="60" v2="62"/>
+ <triangle index="70" v0="61" v1="63" v2="64"/>
+ <triangle index="71" v0="61" v1="62" v2="64"/>
+ <triangle index="72" v0="63" v1="65" v2="66"/>
+ <triangle index="73" v0="63" v1="64" v2="66"/>
+ <triangle index="74" v0="65" v1="67" v2="68"/>
+ <triangle index="75" v0="65" v1="66" v2="68"/>
+ <triangle index="76" v0="32" v1="33" v2="69"/>
+ <triangle index="77" v0="33" v1="69" v2="70"/>
+ <triangle index="78" v0="32" v1="35" v2="71"/>
+ <triangle index="79" v0="32" v1="69" v2="71"/>
+ <triangle index="80" v0="35" v1="37" v2="72"/>
+ <triangle index="81" v0="35" v1="71" v2="72"/>
+ <triangle index="82" v0="37" v1="39" v2="73"/>
+ <triangle index="83" v0="37" v1="72" v2="73"/>
+ <triangle index="84" v0="39" v1="41" v2="74"/>
+ <triangle index="85" v0="39" v1="73" v2="74"/>
+ <triangle index="86" v0="41" v1="43" v2="75"/>
+ <triangle index="87" v0="41" v1="74" v2="75"/>
+ <triangle index="88" v0="43" v1="45" v2="76"/>
+ <triangle index="89" v0="43" v1="75" v2="76"/>
+ <triangle index="90" v0="45" v1="47" v2="77"/>
+ <triangle index="91" v0="45" v1="76" v2="77"/>
+ <triangle index="92" v0="47" v1="48" v2="78"/>
+ <triangle index="93" v0="47" v1="77" v2="78"/>
+ <triangle index="94" v0="48" v1="49" v2="79"/>
+ <triangle index="95" v0="48" v1="78" v2="79"/>
+ <triangle index="96" v0="49" v1="50" v2="80"/>
+ <triangle index="97" v0="49" v1="79" v2="80"/>
+ <triangle index="98" v0="50" v1="52" v2="81"/>
+ <triangle index="99" v0="50" v1="80" v2="81"/>
+ <triangle index="100" v0="52" v1="54" v2="82"/>
+ <triangle index="101" v0="52" v1="81" v2="82"/>
+ <triangle index="102" v0="54" v1="56" v2="83"/>
+ <triangle index="103" v0="54" v1="82" v2="83"/>
+ <triangle index="104" v0="56" v1="58" v2="84"/>
+ <triangle index="105" v0="56" v1="83" v2="84"/>
+ <triangle index="106" v0="58" v1="60" v2="85"/>
+ <triangle index="107" v0="58" v1="84" v2="85"/>
+ <triangle index="108" v0="60" v1="62" v2="86"/>
+ <triangle index="109" v0="60" v1="85" v2="86"/>
+ <triangle index="110" v0="62" v1="64" v2="87"/>
+ <triangle index="111" v0="62" v1="86" v2="87"/>
+ <triangle index="112" v0="64" v1="66" v2="88"/>
+ <triangle index="113" v0="64" v1="87" v2="88"/>
+ <triangle index="114" v0="66" v1="68" v2="89"/>
+ <triangle index="115" v0="66" v1="88" v2="89"/>
+ <triangle index="116" v0="78" v1="79" v2="90"/>
+ <triangle index="117" v0="78" v1="90" v2="91"/>
+ <triangle index="118" v0="79" v1="80" v2="92"/>
+ <triangle index="119" v0="79" v1="90" v2="92"/>
+ <triangle index="120" v0="90" v1="91" v2="93"/>
+ <triangle index="121" v0="91" v1="93" v2="94"/>
+ <triangle index="122" v0="90" v1="92" v2="95"/>
+ <triangle index="123" v0="90" v1="93" v2="95"/>
+ <triangle index="124" v0="93" v1="94" v2="96"/>
+ <triangle index="125" v0="94" v1="96" v2="97"/>
+ <triangle index="126" v0="93" v1="95" v2="98"/>
+ <triangle index="127" v0="93" v1="96" v2="98"/>
+ <triangle index="128" v0="96" v1="97" v2="99"/>
+ <triangle index="129" v0="97" v1="99" v2="100"/>
+ <triangle index="130" v0="96" v1="98" v2="101"/>
+ <triangle index="131" v0="96" v1="99" v2="101"/>
+ <triangle index="132" v0="99" v1="100" v2="102"/>
+ <triangle index="133" v0="100" v1="102" v2="103"/>
+ <triangle index="134" v0="99" v1="101" v2="104"/>
+ <triangle index="135" v0="99" v1="102" v2="104"/>
+ <triangle index="136" v0="102" v1="103" v2="105"/>
+ <triangle index="137" v0="103" v1="105" v2="106"/>
+ <triangle index="138" v0="102" v1="104" v2="107"/>
+ <triangle index="139" v0="102" v1="105" v2="107"/>
+ <triangle index="140" v0="105" v1="106" v2="108"/>
+ <triangle index="141" v0="106" v1="108" v2="109"/>
+ <triangle index="142" v0="105" v1="107" v2="110"/>
+ <triangle index="143" v0="105" v1="108" v2="110"/>
+ <triangle index="144" v0="108" v1="109" v2="111"/>
+ <triangle index="145" v0="109" v1="111" v2="112"/>
+ <triangle index="146" v0="108" v1="110" v2="113"/>
+ <triangle index="147" v0="108" v1="111" v2="113"/>
+ <triangle index="148" v0="111" v1="112" v2="114"/>
+ <triangle index="149" v0="112" v1="114" v2="115"/>
+ <triangle index="150" v0="111" v1="113" v2="116"/>
+ <triangle index="151" v0="111" v1="114" v2="116"/>
+ </cells>
+ <data>
+ <data_entry name="global vertex indices">
+ <meshfunction type="uint" dim="0" size="117">
+ <entity index="0" value="429"/>
+ <entity index="1" value="430"/>
+ <entity index="2" value="471"/>
+ <entity index="3" value="470"/>
+ <entity index="4" value="431"/>
+ <entity index="5" value="472"/>
+ <entity index="6" value="512"/>
+ <entity index="7" value="511"/>
+ <entity index="8" value="513"/>
+ <entity index="9" value="553"/>
+ <entity index="10" value="552"/>
+ <entity index="11" value="554"/>
+ <entity index="12" value="594"/>
+ <entity index="13" value="593"/>
+ <entity index="14" value="595"/>
+ <entity index="15" value="635"/>
+ <entity index="16" value="634"/>
+ <entity index="17" value="636"/>
+ <entity index="18" value="676"/>
+ <entity index="19" value="675"/>
+ <entity index="20" value="677"/>
+ <entity index="21" value="717"/>
+ <entity index="22" value="716"/>
+ <entity index="23" value="718"/>
+ <entity index="24" value="758"/>
+ <entity index="25" value="757"/>
+ <entity index="26" value="759"/>
+ <entity index="27" value="799"/>
+ <entity index="28" value="798"/>
+ <entity index="29" value="800"/>
+ <entity index="30" value="789"/>
+ <entity index="31" value="790"/>
+ <entity index="32" value="831"/>
+ <entity index="33" value="830"/>
+ <entity index="34" value="791"/>
+ <entity index="35" value="832"/>
+ <entity index="36" value="792"/>
+ <entity index="37" value="833"/>
+ <entity index="38" value="793"/>
+ <entity index="39" value="834"/>
+ <entity index="40" value="794"/>
+ <entity index="41" value="835"/>
+ <entity index="42" value="795"/>
+ <entity index="43" value="836"/>
+ <entity index="44" value="796"/>
+ <entity index="45" value="837"/>
+ <entity index="46" value="797"/>
+ <entity index="47" value="838"/>
+ <entity index="48" value="839"/>
+ <entity index="49" value="840"/>
+ <entity index="50" value="841"/>
+ <entity index="51" value="801"/>
+ <entity index="52" value="842"/>
+ <entity index="53" value="802"/>
+ <entity index="54" value="843"/>
+ <entity index="55" value="803"/>
+ <entity index="56" value="844"/>
+ <entity index="57" value="804"/>
+ <entity index="58" value="845"/>
+ <entity index="59" value="805"/>
+ <entity index="60" value="846"/>
+ <entity index="61" value="806"/>
+ <entity index="62" value="847"/>
+ <entity index="63" value="807"/>
+ <entity index="64" value="848"/>
+ <entity index="65" value="808"/>
+ <entity index="66" value="849"/>
+ <entity index="67" value="809"/>
+ <entity index="68" value="850"/>
+ <entity index="69" value="872"/>
+ <entity index="70" value="871"/>
+ <entity index="71" value="873"/>
+ <entity index="72" value="874"/>
+ <entity index="73" value="875"/>
+ <entity index="74" value="876"/>
+ <entity index="75" value="877"/>
+ <entity index="76" value="878"/>
+ <entity index="77" value="879"/>
+ <entity index="78" value="880"/>
+ <entity index="79" value="881"/>
+ <entity index="80" value="882"/>
+ <entity index="81" value="883"/>
+ <entity index="82" value="884"/>
+ <entity index="83" value="885"/>
+ <entity index="84" value="886"/>
+ <entity index="85" value="887"/>
+ <entity index="86" value="888"/>
+ <entity index="87" value="889"/>
+ <entity index="88" value="890"/>
+ <entity index="89" value="891"/>
+ <entity index="90" value="922"/>
+ <entity index="91" value="921"/>
+ <entity index="92" value="923"/>
+ <entity index="93" value="963"/>
+ <entity index="94" value="962"/>
+ <entity index="95" value="964"/>
+ <entity index="96" value="1004"/>
+ <entity index="97" value="1003"/>
+ <entity index="98" value="1005"/>
+ <entity index="99" value="1045"/>
+ <entity index="100" value="1044"/>
+ <entity index="101" value="1046"/>
+ <entity index="102" value="1086"/>
+ <entity index="103" value="1085"/>
+ <entity index="104" value="1087"/>
+ <entity index="105" value="1127"/>
+ <entity index="106" value="1126"/>
+ <entity index="107" value="1128"/>
+ <entity index="108" value="1168"/>
+ <entity index="109" value="1167"/>
+ <entity index="110" value="1169"/>
+ <entity index="111" value="1209"/>
+ <entity index="112" value="1208"/>
+ <entity index="113" value="1210"/>
+ <entity index="114" value="1250"/>
+ <entity index="115" value="1249"/>
+ <entity index="116" value="1251"/>
+ </meshfunction>
+ </data_entry>
+ </data>
+ </mesh>
+</dolfin>
=== added file 'nitsche_method/data/meshes/rotor_cavity.geo'
--- nitsche_method/data/meshes/rotor_cavity.geo 1970-01-01 00:00:00 +0000
+++ nitsche_method/data/meshes/rotor_cavity.geo 2010-03-16 14:39:09 +0000
@@ -0,0 +1,82 @@
+//This file defines start geometry for the driven rotor in the cavity problem.
+
+lb = 0.5;
+lc = 0.3;
+
+// This variable can then be used in the definition of Gmsh's simplest
+// `elementary entity', a `Point'. A Point is defined by a list of
+// four numbers: three coordinates (X, Y and Z), and a characteristic
+// length (lc) that sets the target element size at the point:
+
+//Cavity geometry
+H = 2.0;
+
+Point(1) = {-H, -H, 0, lb};
+Point(2) = {H, -H, 0, lb};
+Point(3) = {H, H, 0, lb};
+Point(4) = {-H, H, 0, lb};
+
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+
+Line Loop(5) = {1,2,3,4};
+
+//Physical Point(1) = {1,2,3,4};
+//Physical Line(1) = {5};
+
+
+//Rotor geometry
+h = 0.1;
+l = 1.0;
+
+//Point(5) = {-h,-h,0,lc};
+//Point(6) = {h,-h,0,lc};
+//Point(7) = {h,h,0,lc};
+//Point(8) = {-h,h,0,lc};
+//
+//Line(6) = {5,6};
+//Line(7) = {6,7};
+//Line(8) = {7,8};
+//Line(9) = {8,5};
+//
+//Line Loop(10) = {6,7,8,9};
+//Plane Surface(101) = {10};
+//Plane Surface(100) = {10};
+//Plane Surface(200) = {5,10};
+
+Point(5) = {h,h,0,lc};
+Point(6) = {h,l,0,lc};
+Point(7) = {-h,l,0,lc};
+Point(8) = {-h,h,0,lc};
+Point(9) = {-l,h,0,lc};
+Point(10) = {-l,-h,0,lc};
+Point(11) = {-h,-h,0,lc};
+Point(12) = {-h,-l,0,lc};
+Point(13) = {h,-l,0,lc};
+Point(14) = {h,-h,0,lc};
+Point(15) = {l,-h,0,lc};
+Point(16) = {l,h,0,lc};
+
+Line(6) = {5,6};
+Line(7) = {6,7};
+Line(8) = {7,8};
+Line(9) = {8,9};
+Line(10) = {9,10};
+Line(11) = {10,11};
+Line(12) = {11,12};
+Line(13) = {12,13};
+Line(14) = {13,14};
+Line(15) = {14,15};
+Line(16) = {14,15};
+Line(17) = {15,16};
+Line(18) = {16,5};
+
+Line Loop(19) = {6,7,8,9,10,11,12,13,14,15,16,17,18};
+
+//Physical Point(2) = {6,7,8,9,10,11,12,13,14,15,16,17,18};
+//Physical Line(2) = {19};
+
+Plane Surface(100) = {19};
+Plane Surface(101) = {5,19};
=== added directory 'nitsche_method/entity_intersection'
=== added file 'nitsche_method/entity_intersection/SConstruct'
--- nitsche_method/entity_intersection/SConstruct 1970-01-01 00:00:00 +0000
+++ nitsche_method/entity_intersection/SConstruct 2010-02-10 03:14:09 +0000
@@ -0,0 +1,16 @@
+# This is an example of a simple SConstruct file for building
+# programs against DOLFIN. To build this demo, just type 'scons'.
+
+import os, commands
+
+# Get compiler from pkg-config
+compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
+
+# Create a SCons Environment based on the main os environment
+env = Environment(ENV=os.environ, CXX=compiler)
+
+# Get compiler flags from pkg-config
+env.ParseConfig('pkg-config --cflags --libs dolfin')
+
+#env.Program(['tetrahedron_triangle_intersection_test.cpp'])
+env.Program(['primitive_intersector_test.cpp'])
=== added file 'nitsche_method/entity_intersection/primitive_intersector_test.cpp'
--- nitsche_method/entity_intersection/primitive_intersector_test.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/entity_intersection/primitive_intersector_test.cpp 2010-03-05 23:36:16 +0000
@@ -0,0 +1,40 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-02-10 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-02-10
+// Last changed: 2010-02-22
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+// =====================================================================================
+
+#include <dolfin/mesh/dolfin_mesh.h>
+#include <iostream>
+
+using namespace dolfin;
+
+int main ()
+{
+ UnitCube cube(3,3,2);
+ cout <<"Total number of cells in Cube:" << cube.num_cells() <<endl;
+
+ UnitSphere sphere(3);
+ cout <<"Total number of cells in Sphere:" << sphere.num_cells() <<endl;
+
+ for (CellIterator cube_cell(cube); !cube_cell.end(); ++cube_cell)
+ {
+ for (CellIterator sphere_cell(cube); !sphere_cell.end(); ++sphere_cell)
+ {
+ if (PrimitiveIntersector::do_intersect(*cube_cell, *sphere_cell))
+ std::cout <<"Cell number " << cube_cell->index() << " intersects with "
+ <<"Cell number " << sphere_cell->index() << std::endl;
+ }
+ }
+
+ return 0;
+}
=== added directory 'nitsche_method/integration_on_complex_domains'
=== added directory 'nitsche_method/integration_on_complex_domains/gauss_projection_green'
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/BaryCenterQuadrature.cpp'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/BaryCenterQuadrature.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/BaryCenterQuadrature.cpp 2010-04-01 10:43:19 +0000
@@ -0,0 +1,270 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-03-17 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-03-17
+// Last changed: 2010-04-01
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//Remark: This is an adapted version of Brian Mirtichs original C code for
+//computin polyhedral mass properties, which can be found at <web-address>
+// =====================================================================================
+
+#include "BaryCenterQuadrature.h"
+#include <math.h>
+
+#define X 0
+#define Y 1
+#define Z 2
+#define SQR(x) ((x)*(x))
+#define CUBE(x) ((x)*(x)*(x))
+
+
+using std::sqrt;
+using namespace dolfin;
+
+
+void BaryCenterQuadrature::compute_quadrature_rule(const Nef_polyhedron_3 & polyhedron)
+{
+ //coordinates of the bary center
+ double bary_center [3] = {0,0,0};
+
+ _weight = 0;
+
+ //determine if polyhedron has a non vanishing 3d measure.
+ //This requires to have volumes marked with one.
+ bool has_volume = false;
+ Nef_polyhedron_3::Volume_const_iterator vi;
+ CGAL_forall_volumes(vi,polyhedron)
+ {
+ std::cout <<"vi->mark: " << vi->mark() <<std::endl;
+ if(vi->mark())
+ {
+ has_volume = true;
+ break;
+ }
+ }
+
+ //variables needing during the computation
+ //components of surface normal
+ double nx, ny, nz;
+ double normal [3] = {0,0,0};
+
+ //start and end point of the edge
+ double source[3] = {0,0,0};
+ double target[3] = {0,0,0};
+
+ //permutation of the coordinates X,Y,Z
+ int A,B,C;
+
+ std::cout <<" Number of volumes :" << polyhedron.number_of_volumes() << std::endl;
+ std::cout << "Outer iteration along faces, next f \n";
+ std::cout <<"++++++++++++++++++++++++++++++++++++++++" << std::endl;
+
+ std::cout << "Polyhedron is 2d Manifold: " << const_cast<Nef_polyhedron_3 &>(polyhedron).is_simple() << std::endl;
+
+ for(Nef_polyhedron_3::Halffacet_const_iterator f = polyhedron.halffacets_begin (),
+ end = polyhedron.halffacets_end();
+ f != end; ++f)
+ {
+ if (f->incident_volume()->mark() || (!has_volume && f->is_twin()))
+ {
+ //compute best projection direction.
+ //get normal of the surface and the offset
+ //FIXME: check right orientation
+ Plane_3 face_plane = f->plane();
+ Vector_3 test_normal = face_plane.orthogonal_vector();
+ nx = CGAL::to_double(face_plane.a());
+ ny = CGAL::to_double(face_plane.b());
+ nz = CGAL::to_double(face_plane.c());
+
+ double normal_length = sqrt(nx*nx + ny*ny + nz*nz);
+ normal[0] = nx/normal_length;
+ normal[1] = ny/normal_length;
+ normal[2] = nz/normal_length;
+
+ nx = fabs(nx);
+ ny = fabs(ny);
+ nz = fabs(nz);
+
+ //Choose permutation to get best conditioned projection direction.
+ if (nx > ny && nx > nz)
+ C = X;
+ else
+ C = (ny > nz) ? Y : Z;
+ A = (C + 1) % 3;
+ B = (A + 1) % 3;
+
+ //begin compute the face integrals
+ //begin compute the projection integrals
+ double a0, a1, da;
+ double b0, b1, db;
+ double a0_2, b0_2;
+ double a0_3, b0_3;
+ double a0_4, b0_4;
+ double a1_2, b1_2;
+ double C1, Ca, Caa, Cb, Cbb, Cab;
+ double Kab;
+ double P1, Pa, Pb, Paa, Pbb, Pab;
+ P1 = Pa = Pb = Paa = Pab = Pbb = 0.0;
+ double Fa, Fb, Fc, Faa, Fbb, Fcc;
+
+ for(Nef_polyhedron_3::Halffacet_cycle_const_iterator fc = f->facet_cycles_begin();
+ fc != f->facet_cycles_end(); ++fc)
+ {
+ //exclude a SHalfloops (whatever this is)
+ if ( fc.is_shalfedge() )
+ {
+ Nef_polyhedron_3::SHalfedge_const_handle h = fc;
+ Nef_polyhedron_3::SHalfedge_around_facet_const_circulator hc(h);
+ Nef_polyhedron_3::SHalfedge_around_facet_const_circulator he(hc);
+ CGAL_For_all(hc,he)
+ {
+ std::cout <<"----------------------------------------" << std::endl;
+
+ Nef_polyhedron_3::SVertex_const_handle v = hc->source();
+
+ const Point_3 & source_point(v->source()->point());
+ const Point_3 & target_point(v->target()->point());
+ std::cout <<"Source point " << source_point << std::endl;
+ std::cout <<"Target point " << target_point <<std::endl;
+
+ source[0] = CGAL::to_double(source_point.x());
+ source[1] = CGAL::to_double(source_point.y());
+ source[2] = CGAL::to_double(source_point.z());
+
+ target[0] = CGAL::to_double(target_point.x());
+ target[1] = CGAL::to_double(target_point.y());
+ target[2] = CGAL::to_double(target_point.z());
+
+ a0 = source[A];
+ b0 = source[B];
+ a1 = target[A];
+ b1 = target[B];
+
+ da = a1 - a0;
+ db = b1 - b0;
+
+ a0_2 = a0 * a0;
+ a0_3 = a0_2 * a0;
+ a0_4 = a0_3 * a0;
+ a1_2 = a1 * a1;
+ b0_2 = b0 * b0;
+ b0_3 = b0_2 * b0;
+ b0_4 = b0_3 * b0;
+ b1_2 = b1 * b1;
+
+ C1 = a1 + a0;
+ Ca = a1*C1 + a0_2;
+ Caa = a1*Ca + a0_3;
+ Cb = b1*(b1 + b0) + b0_2;
+ Cbb = b1*Cb + b0_3;
+ Cab = 3*a1_2 + 2*a1*a0 + a0_2;
+ Kab = a1_2 + 2*a1*a0 + 3*a0_2;
+
+ P1 += db*C1;
+ Pa += db*Ca;
+ Paa += db*Caa;
+ Pb += da*Cb;
+ Pbb += da*Cbb;
+ Pab += db*(b1*Cab + b0*Kab);
+ }
+ std::cout <<"----------------------------------------" << std::endl;
+ }
+ }
+
+ P1 /= 2.0;
+ Pa /= 6.0;
+ Paa /= 12.0;
+ Pb /= -6.0;
+ Pbb /= -12.0;
+ Pab /= 24.0;
+
+ info("P1 : %e",P1);
+ info("Pa : %e",Pa);
+ info("Paa : %e",Paa);
+ info("Pb : %e",Pb);
+ info("Pbb : %e",Pbb);
+ info("Pab : %e",Pab);
+
+ double k1, k2, k3;
+ k1 = 1 / normal[C]; k2 = k1 * k1; k3 = k2 * k1;
+ double w = CGAL::to_double(face_plane.d())/normal_length ;
+ info("Normal[0] : %e ", normal[0]);
+ info("Normal[1] : %e ", normal[1]);
+ info("Normal[2] : %e ", normal[2]);
+ info("w : %e ", w);
+
+ Fa = k1 * Pa;
+ Fb = k1 * Pb;
+ Fc = -k2 * (normal[A]*Pa + normal[B]*Pb + w*P1);
+
+ //debug
+ Faa = k1 * Paa;
+ Fbb = k1 * Pbb;
+ Fcc = k3 * (SQR(normal[A])*Paa + 2*normal[A]*normal[B]*Pab + SQR(normal[B])*Pbb
+ + w*(2*(normal[A]*Pa + normal[B]*Pb) + w*P1));
+
+ info("Fa : %e",Fa);
+ info("Fb : %e",Fb);
+ info("Fc : %e",Fc);
+ info("Faa : %e",Faa);
+ info("Fbb : %e",Fbb);
+ info("Fcc : %e",Fcc);
+ //end computation face integrals
+
+ if (has_volume)
+ {
+ double contribution = (normal[X] * ((A == X) ? Fa : ((B == X) ? Fb : Fc)));
+ info("Contribution: %e", contribution);
+ info("Normal[X]: %e", normal[X]);
+ info("A: %d", A);
+ info("B: %d", B);
+ info("C: %d", C);
+ _weight += normal[X] * ((A == X) ? Fa : ((B == X) ? Fb : Fc));
+ info("Volume: %e", _weight);
+
+ bary_center[A] += normal[A] * Faa;
+ bary_center[B] += normal[B] * Fbb;
+ bary_center[C] += normal[C] * Fcc;
+ }
+ else
+ {
+ info("No volume: P1 is %e",P1);
+ info("No volume: normal[C] is %e", normal[C]);
+ //does not work yet
+ //@todo does not work yet.
+ _weight += -P1/normal[C];
+ info("No volume: weight is %e", _weight);
+ info("Fa : %e",Fa);
+ info("Fb : %e",Fb);
+ info("Fc : %e",Fc);
+ bary_center[A] += -Fa;
+ bary_center[B] += -Fb;
+ bary_center[C] += -Fc;
+ }
+ }
+ }
+
+ //applies also if polyhedron has no no volume, since
+ //we no iterating two times (each halffacet).
+ if (has_volume)
+ {
+ bary_center[X] /= (2*_weight);
+ bary_center[Y] /= (2*_weight);
+ bary_center[Z] /= (2*_weight);
+ }
+ else
+ {
+ bary_center[X] /= _weight;
+ bary_center[Y] /= _weight;
+ bary_center[Z] /= _weight;
+ }
+
+ //use these later directly.
+ _point = Point(3,bary_center);
+}
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/BaryCenterQuadrature.h'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/BaryCenterQuadrature.h 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/BaryCenterQuadrature.h 2010-04-01 10:43:19 +0000
@@ -0,0 +1,69 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-03-17 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-03-17
+// Last changed: 2010-03-30
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+// =====================================================================================
+
+#ifndef __BARYCENTERQUADRATURE_H
+#define __BARYCENTERQUADRATURE_H
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Polyhedron_3.h>
+
+#include <dolfin/mesh/Point.h>
+
+namespace dolfin {
+
+ typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+ typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3;
+ typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron_3;
+ typedef Nef_polyhedron_3::Point_3 Point_3;
+ typedef Nef_polyhedron_3::Plane_3 Plane_3;
+ typedef Nef_polyhedron_3::Vector_3 Vector_3;
+
+ ///This class computes the barycenter of an arbitrary polyhedron or
+ ///polygon and therefore allows for bary center quadrature on complex
+ ///polyhedrons. Note: barycenter quadrature is exact for polynom deg <= 1.
+ class BaryCenterQuadrature {
+
+ public:
+
+ BaryCenterQuadrature() : _weight(0), _point(0,0,0), _length(1) {}
+
+ ///@todo split the ugly implementation up into some subroutines and removes
+ //all the macro etc stuff. Make also working polyhedrons which are
+ //degenerated.
+ void compute_quadrature_rule(const Nef_polyhedron_3 & polyhedron);
+
+ //return pointer
+ ///Gives the points for the last computed Polyhedron.
+ Point * points() {return &_point;}
+ ///Gives the weights for the last computed Polyhedron.
+ double * weights() {return &_weight; }
+
+ ///Number of quadrature points/weights.
+ uint length() {return _length;}
+
+ private:
+
+ double _weight;
+ Point _point;
+ uint _length;
+
+
+ };
+
+}
+
+#endif // ----- #ifndef __BARYCENTERQUADRATURE_H -----
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/README'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/README 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/README 2010-03-18 18:15:26 +0000
@@ -0,0 +1,146 @@
+1. OVERVIEW
+
+ This code accompanies the paper:
+
+ Brian Mirtich, "Fast and Accurate Computation of
+ Polyhedral Mass Properties," journal of graphics
+ tools, volume 1, number 2, 1996.
+
+ It computes the ten volume integrals needed for
+ determining the center of mass, moments of
+ inertia, and products of inertia for a uniform
+ density polyhedron. From this information, a
+ body frame can be computed.
+
+ To compile the program, use an ANSI compiler, and
+ type something like
+
+ % cc volInt.c -O2 -lm -o volInt
+
+
+ Revision history
+
+ 26 Jan 1996 Program creation.
+
+ 3 Aug 1996 Corrected bug arising when polyhedron density
+ is not 1.0. Changes confined to function main().
+ Thanks to Zoran Popovic for catching this one.
+
+
+
+2. POLYHEDRON GEOMETRY FILES
+
+ The program reads a data file specified on the
+ command line. This data file describes the
+ geometry of a polyhedron, and has the following
+ format:
+
+ N
+
+ x_0 y_0 z_0
+ x_1 y_1 z_1
+ .
+ .
+ .
+ x_{N-1} y_{N-1} z_{N-1}
+
+ M
+
+ k1 v_{1,1} v_{1,2} ... v_{1,k1}
+ k2 v_{2,1} v_{2,2} ... v_{2,k2}
+ .
+ .
+ .
+ kM v_{M,1} v_{M,2} ... v_{M,kM}
+
+
+ where:
+ N number of vertices on polyhedron
+ x_i y_i z_i x, y, and z coordinates of ith vertex
+ M number of faces on polyhedron
+ ki number of vertices on ith face
+ v_{i,j} jth vertex on ith face
+
+ In English:
+
+ First the number of vertices are specified. Next
+ the vertices are defined by listing the 3
+ coordinates of each one. Next the number of faces
+ are specified. Finally, the faces themselves are
+ specified. A face is specified by first giving
+ the number of vertices around the polygonal face,
+ followed by the (integer) indices of these
+ vertices. When specifying indices, note that
+ they must be given in counter-clockwise order
+ (when looking at the face from outside the
+ polyhedron), and the vertices are indexed from 0
+ to N-1 for a polyhedron with N faces.
+
+ White space can be inserted (or not) as desired.
+ Three example polyhedron geometry files are included:
+
+ cube A cube, 20 units on a side, centered at
+ the origin and aligned with the coordinate axes.
+
+ tetra A tetrahedron formed by taking the convex
+ hull of the origin, and the points (5,0,0),
+ (0,4,0), and (0,0,3).
+
+ icosa An icosahedron with vertices lying on the unit
+ sphere, centered at the origin.
+
+
+
+3. RUNNING THE PROGRAM
+
+ Simply type,
+
+ % volInt <polyhedron geometry filename>
+
+ The program will read in the geometry of the
+ polyhedron, and the print out the ten volume
+ integrals.
+
+ The program also computes some of the mass
+ properties which may be inferred from the volume
+ integrals. A density of 1.0 is assumed, although
+ this may be changed in function main(). The
+ center of mass is computed as well as the inertia
+ tensor relative to a frame with origin at the
+ center of mass. Note, however, that this matrix
+ may not be diagonal. If not, a diagonalization
+ routine must be performed, to determine the
+ orientation of the true body frame relative to
+ the original model frame. The Jacobi method
+ works quite well (see Numerical Recipes in C by
+ Press, et. al.).
+
+
+
+4. DISCLAIMERS
+
+ 1. The volume integration code has been written
+ to match the development and algorithms presented
+ in the paper, and not with maximum optimization
+ in mind. While inherently very efficient, a few
+ more cycles can be squeezed out of the algorithm.
+ This is left as an exercise. :)
+
+ 2. Don't like global variables? The three
+ procedures which evaluate the volume integrals
+ can be combined into a single procedure with two
+ nested loops. In addition to providing some
+ speedup, all of the global variables can then be
+ made local.
+
+ 3. The polyhedron data structure used by the
+ program is admittedly lame; much better schemes
+ are possible. The idea here is just to give the
+ basic integral evaluation code, which will have
+ to be adjusted for other polyhedron data
+ structures.
+
+ 4. There is no error checking for the input
+ files. Be careful. Note the hard limits
+ #defined for the number of vertices, number of
+ faces, and number of vertices per faces.
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/SConstruct'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/SConstruct 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/SConstruct 2010-04-02 18:46:57 +0000
@@ -0,0 +1,22 @@
+import os, commands
+
+# Import("cppunitCfg")
+
+# Get compiler from pkg-config
+compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
+
+# Create a SCons Environment based on the main os environment
+env = Environment(ENV=os.environ, CXX=compiler)
+
+env.Append(CXXFLAGS = '-ggdb')
+
+# Get compiler flags from pkg-config
+env.ParseConfig('pkg-config --cflags --libs dolfin')
+
+#Program
+env.Program('quadrature_test',['BaryCenterQuadrature.cpp','main.cpp'])
+env.Program('iteration_tetrahedron.cpp')
+env.Program('iteration_tetrahedron_2.cpp')
+
+env.Append(LIBS='-lcppunit')
+env.Program('unittest',['BaryCenterQuadrature.cpp','test.cpp'])
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/alut_snippet.cpp'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/alut_snippet.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/alut_snippet.cpp 2010-03-18 18:15:26 +0000
@@ -0,0 +1,96 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-02-17 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-02-17
+// Last changed: 2010-02-17
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+// =====================================================================================
+
+
+struct
+{
+ Point_3 p1, p2, p3;
+}triangle;
+
+class Shell_explorer
+{
+ bool first;
+ const Nef_polyhedron& N;
+ vector<triangle> m_triangle_list;
+
+public:
+ Shell_explorer(const Nef_polyhedron& N_)
+ : first(true), N(N_) {}
+
+ void visit(Vertex_const_handle v){}
+ void visit(Halfedge_const_handle ){}
+ void visit(Halffacet_const_handle facet)
+ {
+ for (Halffacet_cycle_const_iterator it =
+facet->facet_cycles_begin(); it != facet->facet_cycles_end(); it++)
+ {
+ SHalfedge_const_handle se = SHalfedge_const_handle(it);
+ CGAL_assertion(se!=0);
+ SHalfedge_around_facet_const_circulator hc_start(se);
+ SHalfedge_around_facet_const_circulator hc_end(hc_start);
+ vector <Point_3> vertices;
+ int count = 0;
+ CGAL_For_all(hc_start,hc_end)
+ {
+ SVertex_const_handle svert = hc_start->source();
+ Point_3 vpoint = svert->center_vertex()->point();
+ vertices.push_back(vpoint);
+ }
+ for (int i=0; i<vertices.size()-2; i++)
+ {
+ triangle t;
+ t.p1 = vertices[0];
+ t.p3 = vertices[i+1];
+ t.p2 = vertices[i+2];
+ m_triangle_list.push_back(t);
+ }
+
+ }
+ }
+ void visit(SHalfedge_const_handle ) {}
+ void visit(SHalfloop_const_handle ) {}
+ void visit(SFace_const_handle sf) {}
+
+ vector<triangle>getTriangleList()
+ {
+ vector<triangle> templist;
+ for (int i=0; i<m_triangle_list.size(); i++)
+ templist.push_back(m_triangle_list[i]);
+ m_triangle_list.clear();
+ return templist;
+ }
+ void reset_minimal_vertex() { first = true; }
+};
+
+Usage:
+void main()
+{
+ Nef_polyhedron nPoly;
+ Volume_const_iterator c;
+ Shell_explorer SE(nPoly);
+
+ CGAL_forall_volumes(c,nPoly)
+ {
+ Shell_entry_const_iterator it;
+ int shellcount = 0;
+ vector <triangle> temp;
+ CGAL_forall_shells_of(it,c)
+ {
+ SE.reset_minimal_vertex();
+ entrypaths.visit_shell_objects(SFace_const_handle(it),SE);
+ vector<triangle> tris = SE.getTriangleList();
+ }
+ }
+}
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/cube'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/cube 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/cube 2010-03-18 18:15:26 +0000
@@ -0,0 +1,21 @@
+8
+
+-10 -10 -10
++10 -10 -10
++10 +10 -10
+-10 +10 -10
+-10 -10 +10
++10 -10 +10
++10 +10 +10
+-10 +10 +10
+
+6
+
+4 0 3 2 1
+4 4 5 6 7
+4 0 1 5 4
+4 6 2 3 7
+4 1 2 6 5
+4 0 4 7 3
+
+
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/icosa'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/icosa 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/icosa 2010-03-18 18:15:26 +0000
@@ -0,0 +1,38 @@
+12
+
++0.00000000 +0.00000000 +1.00000000
++0.00000000 +0.00000000 -1.00000000
++0.89442719 +0.00000000 +0.44721360
++0.27639320 +0.85065081 +0.44721360
+-0.72360680 +0.52573111 +0.44721360
+-0.72360680 -0.52573111 +0.44721360
++0.27639320 -0.85065081 +0.44721360
++0.72360680 +0.52573111 -0.44721360
+-0.27639320 +0.85065081 -0.44721360
+-0.89442719 +0.00000000 -0.44721360
+-0.27639320 -0.85065081 -0.44721360
++0.72360680 -0.52573111 -0.44721360
+
+20
+
+3 6 11 2
+3 3 2 7
+3 7 2 11
+3 0 2 3
+3 0 6 2
+3 10 1 11
+3 1 7 11
+3 10 11 6
+3 8 7 1
+3 8 3 7
+3 5 10 6
+3 5 6 0
+3 4 3 8
+3 4 0 3
+3 9 8 1
+3 9 1 10
+3 4 5 0
+3 9 10 5
+3 9 5 4
+3 9 4 8
+
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/iteration_tetrahedron.cpp'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/iteration_tetrahedron.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/iteration_tetrahedron.cpp 2010-03-25 23:33:08 +0000
@@ -0,0 +1,58 @@
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Polyhedron_3.h>
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3;
+typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron_3;
+typedef Nef_polyhedron_3::Point_3 Point_3;
+
+int main()
+{
+
+ Point_3 a(1.0, 0.0, 0.0);
+ Point_3 b(0.0, 1.0, 0.0);
+ Point_3 c(0.0, 0.0, 1.0);
+ Point_3 d(0.0, 0.0, 0.0);
+
+ //Build a tetrahedron
+ Polyhedron_3 P;
+ P.make_tetrahedron(a,b,c,d);
+ Nef_polyhedron_3 NP(P);
+
+ std::cout << "Iteration along face ...\n";
+ int face_counter = 0;
+ for(Nef_polyhedron_3::Halffacet_const_iterator f = NP.halffacets_begin (),
+ end = NP.halffacets_end();
+ f != end;
+ ++f)
+ {
+ if(f->is_twin())
+ continue;
+ std::cout <<"########################################" << std::endl;
+ std::cout << "Face number " << face_counter++ << std::endl;
+ for(Nef_polyhedron_3::Halffacet_cycle_const_iterator fc = f->facet_cycles_begin(),
+ cycles_end = f->facet_cycles_end();
+ fc != cycles_end;
+ ++fc)
+ {
+ if ( fc.is_shalfedge() )
+ {
+ std::cout <<"Edge consist of points: " << std::endl;
+ Nef_polyhedron_3::SHalfedge_const_handle se = fc;
+ Nef_polyhedron_3::SHalfedge_around_facet_const_circulator hc(se);
+ Nef_polyhedron_3::SHalfedge_around_facet_const_circulator he(hc);
+ CGAL_For_all(hc,he)
+ {
+ Nef_polyhedron_3::SVertex_const_handle v = hc->source();
+// Nef_polyhedron_3::Halfedge_const_handle v = hc->source();
+ const Point_3 & source_point(v->source()->point());
+ const Point_3 & target_point(v->target()->point());
+ std::cout <<"Source point " << source_point << std::endl;
+ std::cout <<"Target point " << target_point <<std::endl;
+ }
+ }
+ }
+ }
+ return 0;
+}
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/iteration_tetrahedron_2.cpp'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/iteration_tetrahedron_2.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/iteration_tetrahedron_2.cpp 2010-04-01 10:43:19 +0000
@@ -0,0 +1,140 @@
+#include<vector>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Gmpz.h>
+#include <CGAL/Homogeneous.h>
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3;
+typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron_3;
+ typedef Nef_polyhedron_3::Vector_3 Vector_3;
+typedef Nef_polyhedron_3::Point_3 Point_3;
+typedef Nef_polyhedron_3::Vertex_const_handle Vertex_const_handle;
+typedef Nef_polyhedron_3::Halfedge_const_handle Halfedge_const_handle;
+typedef Nef_polyhedron_3::Halffacet_const_handle Halffacet_const_handle;
+typedef Nef_polyhedron_3::Halffacet_cycle_const_iterator Halffacet_cycle_const_iterator;
+typedef Nef_polyhedron_3::SHalfedge_const_handle SHalfedge_const_handle;
+typedef Nef_polyhedron_3::SHalfedge_around_facet_const_circulator SHalfedge_around_facet_const_circulator;
+typedef Nef_polyhedron_3::SHalfloop_const_handle SHalfloop_const_handle;
+typedef Nef_polyhedron_3::SFace_const_handle SFace_const_handle;
+typedef Nef_polyhedron_3::Volume_const_iterator Volume_const_iterator;
+typedef Nef_polyhedron_3::Shell_entry_const_iterator Shell_entry_const_iterator;
+
+typedef std::vector<Nef_polyhedron_3> PolyhedronList;
+typedef PolyhedronList::const_iterator PolyhedronListIterator;
+
+class Shell_explorer
+{
+ const Nef_polyhedron_3 & N;
+
+public:
+ Shell_explorer(const Nef_polyhedron_3 & N_)
+ : N(N_) {}
+
+ void visit(Vertex_const_handle v) {}
+ void visit(Halfedge_const_handle ){}
+ void visit(Halffacet_const_handle facet)
+ {
+ std::cout << "Visting faces ..." << std::endl;
+ std::cout <<"++++++++++++++++++++++++++++++++++++++++" << std::endl;
+ Vector_3 test_normal = facet->plane().orthogonal_vector();
+ std::cout <<"Normal vector of the plane is " << test_normal <<std::endl;
+
+ for (Halffacet_cycle_const_iterator it =
+ facet->facet_cycles_begin(); it != facet->facet_cycles_end(); it++)
+ {
+ std::cout << "Visting Edges ..." << std::endl;
+ std::cout <<"----------------------------------------" << std::endl;
+ Nef_polyhedron_3::SHalfedge_const_handle h = it;
+ Nef_polyhedron_3::SHalfedge_around_facet_const_circulator hc(h), he(hc);
+ CGAL_For_all(hc,he)
+ {
+ Nef_polyhedron_3::SVertex_const_handle v = hc->source();
+ const Nef_polyhedron_3::Point_3& source = v->source()->point();
+ const Nef_polyhedron_3::Point_3& target = v->target()->point();
+ std::cout <<"Source :" <<source <<std::endl;
+ std::cout <<"Target :" <<target <<std::endl;
+ }
+ }
+ }
+ void visit(SHalfedge_const_handle ) {}
+ void visit(SHalfloop_const_handle ) {}
+ void visit(SFace_const_handle sf) {}
+
+};
+
+int main()
+{
+ Point_3 a1(1.0, 0.0, 0.0);
+ Point_3 b1(0.0, 1.0, 0.0);
+ Point_3 c1(0.0, 0.0, 1.0);
+ Point_3 d1(0.0, 0.0, 0.0);
+
+ Point_3 a2(-1.0, 0.0, 0.0);
+ Point_3 b2(0.0, 1.0, 0.0);
+ Point_3 c2(0.0, 0.0, 1.0);
+ Point_3 d2(0.0, 0.0, 0.0);
+
+ Point_3 a3(1.0, 0.0, 0.0);
+ Point_3 b3(0.0, -1.0, 0.0);
+ Point_3 c3(0.0, 0.0, 1.0);
+ Point_3 d3(0.0, 0.0, 0.0);
+
+ Point_3 a4(-1.0, 0.0, 0.0);
+ Point_3 b4(0.0, -1.0, 0.0);
+ Point_3 c4(0.0, 0.0, 1.0);
+ Point_3 d4(0.0, 0.0, 0.0);
+
+ Point_3 A(2.0, 1.0, 1.0);
+ Point_3 B(1.0, 2.0, 1.0);
+ Point_3 C(1.0, 1.0, 2.0);
+ Point_3 D(1.0, 1.0, 1.0);
+
+ //Build a point tetrahedron
+ Polyhedron_3 P1;
+ P1.make_triangle(a1,b1,c1);
+ Nef_polyhedron_3 N1(P1);
+
+ Polyhedron_3 P2;
+ P2.make_tetrahedron(a2,b2,c2,d2);
+ Nef_polyhedron_3 N2(P2);
+
+ Polyhedron_3 P3;
+ P3.make_tetrahedron(a3,b3,c3,d3);
+ Nef_polyhedron_3 N3(P3);
+
+ Polyhedron_3 P4;
+ P4.make_tetrahedron(a4,b4,c4,d4);
+ Nef_polyhedron_3 N4(P4);
+
+ PolyhedronList polyhedrons;
+
+ polyhedrons.push_back(N1);
+ polyhedrons.push_back(N2);
+
+ //Build a point tetrahedron
+ for (PolyhedronListIterator pi = polyhedrons.begin(), end = polyhedrons.end();
+ pi != end; ++pi)
+ {
+ Volume_const_iterator vi;
+ Shell_explorer SE(*pi);
+
+ std::cout <<"Next polyhedron:\n"
+ << "##############################" <<std::endl;
+ int ic = 0;
+ CGAL_forall_volumes(vi,*pi)
+ {
+ std::cout << "Volume " << ic++ << " has mark " << vi->mark()<< std::endl;
+
+ Shell_entry_const_iterator it;
+ CGAL_forall_shells_of(it,vi)
+ {
+ pi->visit_shell_objects(SFace_const_handle(it),SE);
+ }
+ std::cout << "\n\n";
+ }
+ }
+ return 0;
+}
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/main.cpp'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/main.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/main.cpp 2010-04-01 10:43:19 +0000
@@ -0,0 +1,190 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-03-19 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-03-19
+// Last changed: 2010-03-31
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+// =====================================================================================
+
+#include "BaryCenterQuadrature.h"
+#include <dolfin.h>
+
+#include <iostream>
+#include<vector>
+
+using namespace dolfin;
+
+typedef std::vector<Nef_polyhedron_3> PolyhedronList;
+typedef PolyhedronList::const_iterator PolyhedronListIterator;
+
+
+using namespace dolfin;
+
+int main()
+{
+
+ Point_3 e0(0.0, 0.0, 0.0);
+
+ Point_3 e1(1.0, 0.0, 0.0);
+ Point_3 e2(0.0, 1.0, 0.0);
+ Point_3 e3(0.0, 0.0, 1.0);
+
+ Point_3 _e1(-1.0, 0.0, 0.0);
+ Point_3 _e2(0.0, -1.0, 0.0);
+ Point_3 _e3(0.0, 0.0, -1.0);
+
+ Point_3 A(2.0, 1.0, 1.0);
+ Point_3 B(1.0, 2.0, 1.0);
+ Point_3 C(1.0, 1.0, 2.0);
+ Point_3 D(1.0, 1.0, 1.0);
+
+ //Create tetrahedrons with unit vectors
+ Polyhedron_3 P1;
+ P1.make_tetrahedron(e0,e1,e2,e3);
+ Nef_polyhedron_3 N1(P1);
+
+ Polyhedron_3 P2;
+ P2.make_tetrahedron(e0,_e1,e2,e3);
+ Nef_polyhedron_3 N2(P2);
+
+ Polyhedron_3 P3;
+ P3.make_tetrahedron(e0,e1,_e2,e3);
+ Nef_polyhedron_3 N3(P3);
+
+ Polyhedron_3 P4;
+ P4.make_tetrahedron(e0,_e1,_e2,e3);
+ Nef_polyhedron_3 N4(P4);
+
+ Polyhedron_3 P5;
+ P5.make_tetrahedron(e0,e1,e2,_e3);
+ Nef_polyhedron_3 N5(P5);
+
+ Polyhedron_3 P6;
+ P6.make_tetrahedron(e0,_e1,e2,_e3);
+ Nef_polyhedron_3 N6(P6);
+
+ Polyhedron_3 P7;
+ P7.make_tetrahedron(e0,e1,_e2,_e3);
+ Nef_polyhedron_3 N7(P7);
+
+ Polyhedron_3 P8;
+ P8.make_tetrahedron(e0,_e1,_e2,_e3);
+ Nef_polyhedron_3 N8(P8);
+
+ //Create triangles
+ Polyhedron_3 P9;
+ P9.make_triangle(e0,e1,e2);
+ Nef_polyhedron_3 N9(P9);
+
+ Polyhedron_3 P10;
+ P10.make_triangle(e0,e1,e3);
+ Nef_polyhedron_3 N10(P10);
+
+ Polyhedron_3 P11;
+ P11.make_triangle(e0,e2,e3);
+ Nef_polyhedron_3 N11(P11);
+
+ Polyhedron_3 P12;
+ P12.make_triangle(e1,e2,e3);
+ Nef_polyhedron_3 N12(P12);
+
+ Polyhedron_3 P13;
+ P13.make_triangle(e0,_e1,e2);
+ Nef_polyhedron_3 N13(P13);
+
+ Polyhedron_3 P14;
+ P14.make_triangle(e0,_e1,e3);
+ Nef_polyhedron_3 N14(P14);
+
+ Polyhedron_3 P15;
+ P15.make_triangle(e0,e2,e3);
+ Nef_polyhedron_3 N15(P15);
+
+ Polyhedron_3 P16;
+ P16.make_triangle(_e1,e2,e3);
+ Nef_polyhedron_3 N16(P16);
+
+ PolyhedronList polyhedrons;
+ PolyhedronList polygons;
+
+ //single polyhedrons
+ polyhedrons.push_back(N1);
+ polyhedrons.push_back(N2);
+ polyhedrons.push_back(N3);
+ polyhedrons.push_back(N4);
+ polyhedrons.push_back(N5);
+ polyhedrons.push_back(N6);
+ polyhedrons.push_back(N7);
+ polyhedrons.push_back(N8);
+
+ //combined tetrehedrons
+ polyhedrons.push_back(N1 + N2);
+ polyhedrons.push_back(N1 + N2 + N3);
+ polyhedrons.push_back(N1 + N2 + N3 + N4);
+
+ polyhedrons.push_back(N5 + N6);
+ polyhedrons.push_back(N5 + N6 + N7);
+ polyhedrons.push_back(N5 + N6 + N7 + N9);
+
+ polyhedrons.push_back(N1 + N2 + N3 + N4 + N5 + N6 + N7 + N4);
+
+ //polygons
+ polygons.push_back(N9);
+ polygons.push_back(N10);
+ polygons.push_back(N11);
+ polygons.push_back(N12);
+ polygons.push_back(N13);
+ polygons.push_back(N14);
+ polygons.push_back(N15);
+ polygons.push_back(N16);
+
+ //combined polygons
+ polygons.push_back(N9 + N10);
+ polygons.push_back(N9 + N10 + N11);
+ polygons.push_back(N9 + N10 + N11 + N12);
+
+ polygons.push_back(N13 + N14);
+ polygons.push_back(N13 + N14 + N15);
+ polygons.push_back(N13 + N14 + N15 + N16);
+
+ polygons.push_back(N9 + N10 + N11 + N12 + N13 + N14 + N15 + N12);
+
+ BaryCenterQuadrature quadrature_rule;
+
+ int counter = 0;
+ std::cout << "Quadrature rule for Polyhdrons\n" ;
+ std::cout <<"------------------------------\n";
+ for (PolyhedronListIterator pi = polyhedrons.begin(), end = polyhedrons.end();
+ pi != end; ++pi)
+ {
+ cout <<"Next polyhedron: " << ++counter << endl
+ << "##############################" <<endl;
+ quadrature_rule.compute_quadrature_rule(*pi);
+ cout <<"Volume: " << *(quadrature_rule.weights()) <<endl;
+ cout <<"Barycenter: "
+ << quadrature_rule.points()[0] << endl <<endl;
+ }
+
+ counter = 0;
+ std::cout << "Quadrature rule for Polygons\n" ;
+ std::cout <<"------------------------------\n";
+ for (PolyhedronListIterator pi = polygons.begin(), end = polygons.end();
+ pi != end; ++pi)
+ {
+ cout <<"Next polygons: " << ++counter << endl
+ << "##############################" <<endl;
+ quadrature_rule.compute_quadrature_rule(*pi);
+ cout <<"Volume: " << *(quadrature_rule.weights()) <<endl;
+ cout <<"Barycenter: "
+ << quadrature_rule.points()[0] << endl <<endl;
+ }
+
+ return 0;
+}
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/martin_baeker_snippet.cpp'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/martin_baeker_snippet.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/martin_baeker_snippet.cpp 2010-03-18 18:15:26 +0000
@@ -0,0 +1,93 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-02-17 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-02-17
+// Last changed: 2010-02-17
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+// =====================================================================================
+
+ for(Nef_polyhedron::Halffacet_const_iterator
+ f = Cube.halffacets_begin (),
+ end = Cube.halffacets_end();
+ f != end; ++f)
+ {
+ if(f->is_twin()) continue;
+ std::cout << "Outer iteration, next f \n";
+ for(Nef_polyhedron::Halffacet_cycle_const_iterator
+ fc = f->facet_cycles_begin(),
+ end = f->facet_cycles_end();
+ fc != end; ++fc)
+ {
+ std::cout << "Inner iteration, next fc \n";
+ if ( fc.is_shalfedge() )
+ {
+ std::cout << "Halfedge consists of points \n";
+
+ Nef_polyhedron::SHalfedge_const_handle h = fc;
+ Nef_polyhedron::SHalfedge_around_facet_const_circulator hc(h), he(hc);
+ CGAL_For_all(hc,he)
+ { // all vertex coordinates in facet cycle
+ Nef_polyhedron::SVertex_const_handle v = hc->source();
+ const Nef_polyhedron::Point_3& point = v->source()->point();
+ std::cout << "p: " << CGAL::to_double(point.x()) << " " <<
+ CGAL::to_double(point.y()) << " " <<
+ CGAL::to_double(point.z()) << std::endl;
+
+ }
+
+ }
+ }
+
+ }
+
+ // Alternative to iterating over halfedges:
+// Nef_polyhedron::Halfedge_const_iterator e = Cube.halfedges_begin();
+// CGAL_forall_halfedges(e,Cube) {
+
+
+ for(Nef_polyhedron::Halfedge_const_iterator
+ e = Cube.halfedges_begin(),
+ end = Cube.halfedges_end();
+ e != end; ++e)
+ {
+ const Nef_polyhedron::Vertex_const_handle& s = e->source();
+ const Nef_polyhedron::Vertex_const_handle& t = e->twin()->source();
+ const Nef_polyhedron::Point_3& a = s->point();
+ const Nef_polyhedron::Point_3& b = t->point();
+ std::cout << "From a: " << CGAL::to_double(a.x()) << " " <<
+ CGAL::to_double(a.y()) << " " <<
+ CGAL::to_double(a.z()) << " to b: " <<
+ CGAL::to_double(b.x()) << " " <<
+ CGAL::to_double(b.y()) << " " <<
+ CGAL::to_double(b.z()) <<std::endl;
+ typedef Nef_polyhedron::Object_handle Object_handle;
+ typedef Nef_polyhedron::Vertex_const_handle Vertex_const_handle;
+ Vertex_const_handle v;
+ Object_handle o = Cube.locate(a);
+ if(CGAL::assign(v,o))
+ std::cout << " a is a vertex" << std::endl;
+ else
+ std::cout << " a is not a vertex" << std::endl;
+ }
+
+// Nef_polyhedron::Vertex_const_iterator v;
+
+// CGAL_forall_vertices(v, Cube.sncp())
+ for(Nef_polyhedron::Vertex_const_iterator
+ v = Cube.vertices_begin(),
+ end = Cube.vertices_end();
+ v != end; ++v)
+ {
+ const Nef_polyhedron::Point_3& a = v->point();
+ std::cout << "Vertex: " << CGAL::to_double(a.x()) << " " <<
+ CGAL::to_double(a.y()) << " " <<
+ CGAL::to_double(a.z()) << std::endl;
+
+ }
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/test.cpp'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/test.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/test.cpp 2010-04-02 22:47:40 +0000
@@ -0,0 +1,598 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-04-01 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-04-01
+// Last changed: 2010-04-03
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+//Description: Unittest for BaryCenterQuadrature. =====================================================================================
+
+#include<vector>
+#include <iostream>
+
+#include <dolfin.h>
+#include <dolfin/common/unittest.h>
+
+#include <CGAL/Nef_polyhedron_3.h>
+#include "BaryCenterQuadrature.h"
+
+using namespace dolfin;
+
+typedef Nef_polyhedron_3::Aff_transformation_3 Aff_transformation_3;
+
+typedef std::vector<Nef_polyhedron_3> PolyhedronList;
+typedef PolyhedronList::const_iterator PolyhedronListIterator;
+
+typedef std::vector<int> IntList;
+typedef std::vector<int>::const_iterator IntListIterator;
+typedef std::vector<double> DoubleList;
+typedef std::vector<double>::const_iterator DoubleListIterator;
+typedef std::vector<Point> PointList;
+typedef std::vector<Point>::const_iterator PointListIterator;
+
+class BaryCenter : public CppUnit::TestFixture
+{
+ CPPUNIT_TEST_SUITE(BaryCenter);
+ CPPUNIT_TEST(testSimplePolyhedrons);
+ CPPUNIT_TEST(testSimplePolygons);
+ CPPUNIT_TEST(testComplexPolyhedrons);
+ CPPUNIT_TEST(testComplexPolygons);
+ CPPUNIT_TEST_SUITE_END();
+
+
+ void almost_equal_points(const Point & p1, const Point & p2, double delta)
+ {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(p1.x(),p2.x(),delta);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(p1.y(),p2.y(),delta);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(p1.z(),p2.z(),delta);
+ }
+
+ //Helper function to create reference polyhedrons.
+ void add_test_polyhedron(const Point_3 & p1,
+ const Point_3 & p2,
+ const Point_3 & p3,
+ const Point_3 & p4,
+ PolyhedronList & polyhedrons
+ )
+ {
+ Polyhedron_3 P;
+ P.make_tetrahedron(p1,p2,p3,p4);
+ Nef_polyhedron_3 N(P);
+ polyhedrons.push_back(N);
+ }
+
+ void add_test_polyhedron(const Point_3 & p1,
+ const Point_3 & p2,
+ const Point_3 & p3,
+ PolyhedronList & polyhedrons
+ )
+ {
+ Polyhedron_3 P;
+ P.make_triangle(p1,p2,p3);
+ Nef_polyhedron_3 N(P);
+ polyhedrons.push_back(N);
+ }
+
+ //Helper function to union disjoint polyhedrons and to compute the volume and barycenter.
+ //Indices indicate which polyhedrons should be unioned. No checks at all (index, disjointness etc)
+ //Computed polyhedrons, volumes and barycenters will be append to the given list.
+ void add_disjoint_polyhedrons(const IntList indices,
+ PolyhedronList & polyhedrons,
+ DoubleList & volumes,
+ PointList & points)
+ {
+ double volume = 0;
+ Point point(0,0,0);
+ Nef_polyhedron_3 polyhedron;
+
+ for (IntListIterator i = indices.begin(); i != indices.end(); ++i)
+ {
+ polyhedron += polyhedrons[*i];
+ point += volumes[*i] * points[*i] ;
+ volume += volumes[*i];
+ }
+ point /= volume;
+
+ polyhedrons.push_back(polyhedron);
+ volumes.push_back(volume);
+ points.push_back(point);
+ }
+
+
+ void testSimplePolyhedrons()
+ {
+ //Create origin and unit vectors.
+ Point_3 e0(0.0, 0.0, 0.0);
+
+ Point_3 e1(1.0, 0.0, 0.0);
+ Point_3 e2(0.0, 1.0, 0.0);
+ Point_3 e3(0.0, 0.0, 1.0);
+
+ Point_3 _e1(-1.0, 0.0, 0.0);
+ Point_3 _e2(0.0, -1.0, 0.0);
+ Point_3 _e3(0.0, 0.0, -1.0);
+
+ //Create tetrahedrons with unit vectors and reference results.
+ DoubleList reference_volumes;
+ PointList reference_bary_centers;
+ PolyhedronList reference_polyhedrons;
+
+ add_test_polyhedron(e0,e1,e2,e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(0.25,0.25,0.25));
+
+ add_test_polyhedron(e0,_e1,e2,e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(-0.25,0.25,0.25));
+
+ add_test_polyhedron(e0,e1,_e2,e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(0.25,-0.25,0.25));
+
+ add_test_polyhedron(e0,_e1,_e2,e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(-0.25,-0.25,0.25));
+
+ add_test_polyhedron(e0,e1,e2,_e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(0.25,0.25,-0.25));
+
+ add_test_polyhedron(e0,_e1,e2,_e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(-0.25,0.25,-0.25));
+
+ add_test_polyhedron(e0,e1,_e2,_e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(0.25,-0.25,-0.25));
+
+ add_test_polyhedron(e0,_e1,_e2,_e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(-0.25,-0.25,-0.25));
+
+ //Add sum of polyhedrons
+ IntList add_indices;
+ add_indices.push_back(0);
+
+ add_indices.push_back(1);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(2);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(3);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(4);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(5);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(6);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(7);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ //Add translated version
+ //Upper halfspace
+ Nef_polyhedron_3 polyhedron = reference_polyhedrons[0];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(1, 1, 1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(1.25,1.25,1.25));
+
+ polyhedron = reference_polyhedrons[1];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(-1, 1, 1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(-1.25,1.25,1.25));
+
+ polyhedron = reference_polyhedrons[2];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(1, -1, 1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(1.25,-1.25,1.25));
+
+ polyhedron = reference_polyhedrons[3];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(-1, -1, 1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(-1.25,-1.25,1.25));
+
+ polyhedron = reference_polyhedrons[4];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(1, 1, -1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(1.25,1.25,-1.25));
+
+ polyhedron = reference_polyhedrons[5];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(-1, 1, -1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(-1.25,1.25,-1.25));
+
+ polyhedron = reference_polyhedrons[6];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(1, -1, -1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(1.25,-1.25,-1.25));
+
+ polyhedron = reference_polyhedrons[7];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(-1, -1, -1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(1.0/6.0);
+ reference_bary_centers.push_back(Point(-1.25,-1.25,-1.25));
+
+ //Add the disjoint union of the translated polyhedrons.
+ add_indices.clear();
+ add_indices.push_back(15);
+
+ add_indices.push_back(16);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(17);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(18);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(19);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(20);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(21);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(22);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ //Instantiate quadrature rule
+ BaryCenterQuadrature quadrature_rule;
+
+ //Check volume and barycenter for polyhedrons
+ for (dolfin::uint i = 0; i < reference_polyhedrons.size(); ++i)
+ {
+ quadrature_rule.compute_quadrature_rule(reference_polyhedrons[i]);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(reference_volumes[i],
+ quadrature_rule.weights()[0], 1.0e-12);
+ almost_equal_points(reference_bary_centers[i],
+ quadrature_rule.points()[0], 1.0e-12);
+ }
+ }
+
+ void testSimplePolygons()
+ {
+ //Create origin and unit vectors.
+ Point_3 e0(0.0, 0.0, 0.0);
+
+ Point_3 e1(1.0, 0.0, 0.0);
+ Point_3 e2(0.0, 1.0, 0.0);
+ Point_3 e3(0.0, 0.0, 1.0);
+
+ Point_3 _e1(-1.0, 0.0, 0.0);
+ Point_3 _e2(0.0, -1.0, 0.0);
+ Point_3 _e3(0.0, 0.0, -1.0);
+
+ //Create tetrahedrons with unit vectors and reference results.
+ DoubleList reference_volumes;
+ PointList reference_bary_centers;
+ PolyhedronList reference_polyhedrons;
+
+ //skew plane upper e3 plane
+ add_test_polyhedron(e1,e2,e3,reference_polyhedrons);
+ //todo find exact values
+ reference_volumes.push_back(8.660254e-01);
+ reference_bary_centers.push_back(Point(1.0/3.0,1.0/3.0,1.0/3.0));
+
+ add_test_polyhedron(_e1,e2,e3,reference_polyhedrons);
+ //todo find exact values
+ reference_volumes.push_back(8.660254e-01);
+ reference_bary_centers.push_back(Point(-1.0/3.0,1.0/3.0,1.0/3.0));
+
+ add_test_polyhedron(_e1,_e2,e3,reference_polyhedrons);
+ //todo find exact values
+ reference_volumes.push_back(8.660254e-01);
+ reference_bary_centers.push_back(Point(-1.0/3.0,-1.0/3.0,1.0/3.0));
+
+ add_test_polyhedron(e1,_e2,e3,reference_polyhedrons);
+ //todo find exact values
+ reference_volumes.push_back(8.660254e-01);
+ reference_bary_centers.push_back(Point(1.0/3.0,-1.0/3.0,1.0/3.0));
+
+ //skew plane lower -e3 plane
+ add_test_polyhedron(e1,e2,_e3,reference_polyhedrons);
+ //todo find exact values
+ reference_volumes.push_back(8.660254e-01);
+ reference_bary_centers.push_back(Point(1.0/3.0,1.0/3.0,-1.0/3.0));
+
+ add_test_polyhedron(_e1,e2,_e3,reference_polyhedrons);
+ //todo find exact values
+ reference_volumes.push_back(8.660254e-01);
+ reference_bary_centers.push_back(Point(-1.0/3.0,1.0/3.0,-1.0/3.0));
+
+ add_test_polyhedron(_e1,_e2,_e3,reference_polyhedrons);
+ //todo find exact values
+ reference_volumes.push_back(8.660254e-01);
+ reference_bary_centers.push_back(Point(-1.0/3.0,-1.0/3.0,-1.0/3.0));
+
+ add_test_polyhedron(e1,_e2,_e3,reference_polyhedrons);
+ //todo find exact values
+ reference_volumes.push_back(8.660254e-01);
+ reference_bary_centers.push_back(Point(1.0/3.0,-1.0/3.0,-1.0/3.0));
+
+ //e1-e2 plane
+ add_test_polyhedron(e0,e1,e2,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(1.0/3.0,1.0/3.0,0.0));
+
+ add_test_polyhedron(e0,_e1,e2,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(-1.0/3.0,1.0/3.0,0.0));
+
+ add_test_polyhedron(e0,_e1,_e2,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(-1.0/3.0,-1.0/3.0,0.0));
+
+ add_test_polyhedron(e0,e1,_e2,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(1.0/3.0,-1.0/3.0,0.0));
+
+ //e1-e3 plane
+ add_test_polyhedron(e0,e1,e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(1.0/3.0,0.0,1.0/3.0));
+
+ add_test_polyhedron(e0,_e1,e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(-1.0/3.0,0.0,1.0/3.0));
+
+ add_test_polyhedron(e0,_e1,_e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(-1.0/3.0,0.0,-1.0/3.0));
+
+ add_test_polyhedron(e0,e1,_e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(1.0/3.0,0.0,-1.0/3.0));
+
+ //e2-e3 plane
+ add_test_polyhedron(e0,e2,e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(0.0,1.0/3.0,1.0/3.0));
+
+ add_test_polyhedron(e0,_e2,e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(0.0,-1.0/3.0,1.0/3.0));
+
+ add_test_polyhedron(e0,_e2,_e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(0.0,-1.0/3.0,-1.0/3.0));
+
+ add_test_polyhedron(e0,e2,_e3,reference_polyhedrons);
+ reference_volumes.push_back(1.0/2.0);
+ reference_bary_centers.push_back(Point(0.0,1.0/3.0,-1.0/3.0));
+
+ //Test sum of polyhedrons
+ IntList add_indices;
+ add_indices.push_back(0);
+
+
+ add_indices.push_back(1);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(2);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(3);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(4);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(5);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(6);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(7);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(8);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(9);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+
+ add_indices.push_back(10);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(11);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(12);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(13);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(14);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ add_indices.push_back(15);
+ add_disjoint_polyhedrons(add_indices,
+ reference_polyhedrons,
+ reference_volumes,
+ reference_bary_centers);
+
+ Nef_polyhedron_3 polyhedron = reference_polyhedrons[0];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(1, 1, 1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(reference_volumes[0]);
+ reference_bary_centers.push_back(reference_bary_centers[0] + Point(1,1,1));
+
+ polyhedron = reference_polyhedrons[1];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(-1, 1, 1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(reference_volumes[1]);
+ reference_bary_centers.push_back(reference_bary_centers[1] + Point(-1,1,1));
+
+ polyhedron = reference_polyhedrons[2];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(1, -1, 1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(reference_volumes[2]);
+ reference_bary_centers.push_back(reference_bary_centers[2] + Point(1,-1,1));
+
+ polyhedron = reference_polyhedrons[3];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(-1, -1, 1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(reference_volumes[3]);
+ reference_bary_centers.push_back(reference_bary_centers[3] + Point(-1,-1,1));
+
+ polyhedron = reference_polyhedrons[4];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(1, 1, -1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(reference_volumes[4]);
+ reference_bary_centers.push_back(reference_bary_centers[4] + Point(1,1,-1));
+
+ polyhedron = reference_polyhedrons[5];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(-1, 1, -1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(reference_volumes[5]);
+ reference_bary_centers.push_back(reference_bary_centers[5] + Point(-1,1,-1));
+
+ polyhedron = reference_polyhedrons[6];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(1, -1, -1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(reference_volumes[6]);
+ reference_bary_centers.push_back(reference_bary_centers[6] + Point(1,-1,-1));
+
+ polyhedron = reference_polyhedrons[7];
+ polyhedron.transform(Aff_transformation_3(CGAL::TRANSLATION, Vector_3(-1, -1, -1)));
+ reference_polyhedrons.push_back(polyhedron);
+ reference_volumes.push_back(reference_volumes[7]);
+ reference_bary_centers.push_back(reference_bary_centers[7] + Point(-1,-1,-1));
+
+ //Instantiate quadrature rule
+ BaryCenterQuadrature quadrature_rule;
+
+ //Check volume and barycenter for polyhedrons
+ for (dolfin::uint i = 0; i < reference_polyhedrons.size(); ++i)
+ {
+ quadrature_rule.compute_quadrature_rule(reference_polyhedrons[i]);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(reference_volumes[i],
+ quadrature_rule.weights()[0], 1.0e-5);
+ almost_equal_points(reference_bary_centers[i],
+ quadrature_rule.points()[0], 1.0e-5);
+ }
+ }
+
+ void testComplexPolyhedrons()
+ {
+ }
+
+ void testComplexPolygons()
+ {
+ }
+
+};
+
+CPPUNIT_TEST_SUITE_REGISTRATION(BaryCenter);
+
+int main()
+{
+ DOLFIN_TEST;
+}
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/tetra'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/tetra 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/tetra 2010-03-18 18:15:26 +0000
@@ -0,0 +1,16 @@
+4
+
+0 0 0
+5 0 0
+0 4 0
+0 0 3
+
+4
+
+3 0 3 2
+3 3 0 1
+3 2 1 0
+3 1 2 3
+
+
+
=== added file 'nitsche_method/integration_on_complex_domains/gauss_projection_green/volInt.c'
--- nitsche_method/integration_on_complex_domains/gauss_projection_green/volInt.c 1970-01-01 00:00:00 +0000
+++ nitsche_method/integration_on_complex_domains/gauss_projection_green/volInt.c 2010-03-26 13:05:51 +0000
@@ -0,0 +1,383 @@
+
+
+ /*******************************************************
+ * *
+ * volInt.c *
+ * *
+ * This code computes volume integrals needed for *
+ * determining mass properties of polyhedral bodies. *
+ * *
+ * For more information, see the accompanying README *
+ * file, and the paper *
+ * *
+ * Brian Mirtich, "Fast and Accurate Computation of *
+ * Polyhedral Mass Properties," journal of graphics *
+ * tools, volume 1, number 1, 1996. *
+ * *
+ * This source code is public domain, and may be used *
+ * in any way, shape or form, free of charge. *
+ * *
+ * Copyright 1995 by Brian Mirtich *
+ * *
+ * mirtich@xxxxxxxxxxxxxxx *
+ * http://www.cs.berkeley.edu/~mirtich *
+ * *
+ *******************************************************/
+
+/*
+ Revision history
+
+ 26 Jan 1996 Program creation.
+
+ 3 Aug 1996 Corrected bug arising when polyhedron density
+ is not 1.0. Changes confined to function main().
+ Thanks to Zoran Popovic for catching this one.
+
+ 27 May 1997 Corrected sign error in translation of inertia
+ product terms to center of mass frame. Changes
+ confined to function main(). Thanks to
+ Chris Hecker.
+*/
+
+
+
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+using namespace std;
+
+/*
+ ============================================================================
+ constants
+ ============================================================================
+*/
+
+#define MAX_VERTS 100 /* maximum number of polyhedral vertices */
+#define MAX_FACES 100 /* maximum number of polyhedral faces */
+#define MAX_POLYGON_SZ 10 /* maximum number of verts per polygonal face */
+
+#define X 0
+#define Y 1
+#define Z 2
+
+/*
+ ============================================================================
+ macros
+ ============================================================================
+*/
+
+#define SQR(x) ((x)*(x))
+#define CUBE(x) ((x)*(x)*(x))
+
+/*
+ ============================================================================
+ data structures
+ ============================================================================
+*/
+
+typedef struct {
+ int numVerts;
+ double norm[3];
+ double w;
+ int verts[MAX_POLYGON_SZ];
+ struct polyhedron *poly;
+} FACE;
+
+typedef struct polyhedron {
+ int numVerts, numFaces;
+ double verts[MAX_VERTS][3];
+ FACE faces[MAX_FACES];
+} POLYHEDRON;
+
+
+/*
+ ============================================================================
+ globals
+ ============================================================================
+*/
+
+static int A; /* alpha */
+static int B; /* beta */
+static int C; /* gamma */
+
+/* projection integrals */
+static double P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb;
+
+/* face integrals */
+static double Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca;
+
+/* volume integrals */
+static double T0, T1[3], T2[3], TP[3];
+
+
+/*
+ ============================================================================
+ read in a polyhedron
+ ============================================================================
+*/
+
+void readPolyhedron(char *name, POLYHEDRON *p)
+{
+ FILE *fp;
+ char line[200], *c;
+ int i, j, n;
+ double dx1, dy1, dz1, dx2, dy2, dz2, nx, ny, nz, len;
+ FACE *f;
+
+
+/* if (!(fp = fopen(name, "r"))) {*/
+/* printf("i/o error\n");*/
+/* exit(1);*/
+ }
+
+ fscanf(fp, "%d", &p->numVerts);
+ printf("Reading in %d vertices\n", p->numVerts);
+ for (i = 0; i < p->numVerts; i++)
+ fscanf(fp, "%lf %lf %lf",
+ &p->verts[i][X], &p->verts[i][Y], &p->verts[i][Z]);
+
+ fscanf(fp, "%d", &p->numFaces);
+ printf("Reading in %d faces\n", p->numFaces);
+ for (i = 0; i < p->numFaces; i++) {
+ f = &p->faces[i];
+ f->poly = p;
+ fscanf(fp, "%d", &f->numVerts);
+ for (j = 0; j < f->numVerts; j++) fscanf(fp, "%d", &f->verts[j]);
+
+ /* compute face normal and offset w from first 3 vertices */
+ dx1 = p->verts[f->verts[1]][X] - p->verts[f->verts[0]][X];
+ dy1 = p->verts[f->verts[1]][Y] - p->verts[f->verts[0]][Y];
+ dz1 = p->verts[f->verts[1]][Z] - p->verts[f->verts[0]][Z];
+ dx2 = p->verts[f->verts[2]][X] - p->verts[f->verts[1]][X];
+ dy2 = p->verts[f->verts[2]][Y] - p->verts[f->verts[1]][Y];
+ dz2 = p->verts[f->verts[2]][Z] - p->verts[f->verts[1]][Z];
+ nx = dy1 * dz2 - dy2 * dz1;
+ ny = dz1 * dx2 - dz2 * dx1;
+ nz = dx1 * dy2 - dx2 * dy1;
+ len = sqrt(nx * nx + ny * ny + nz * nz);
+ f->norm[X] = nx / len;
+ f->norm[Y] = ny / len;
+ f->norm[Z] = nz / len;
+ f->w = - f->norm[X] * p->verts[f->verts[0]][X]
+ - f->norm[Y] * p->verts[f->verts[0]][Y]
+ - f->norm[Z] * p->verts[f->verts[0]][Z];
+
+ }
+
+ fclose(fp);
+
+}
+
+/*
+ ============================================================================
+ compute mass properties
+ ============================================================================
+*/
+
+
+/* compute various integrations over projection of face */
+void compProjectionIntegrals(FACE *f)
+{
+ double a0, a1, da;
+ double b0, b1, db;
+ double a0_2, a0_3, a0_4, b0_2, b0_3, b0_4;
+ double a1_2, a1_3, b1_2, b1_3;
+ double C1, Ca, Caa, Caaa, Cb, Cbb, Cbbb;
+ double Cab, Kab, Caab, Kaab, Cabb, Kabb;
+ int i;
+
+ P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0;
+
+ for (i = 0; i < f->numVerts; i++) {
+ a0 = f->poly->verts[f->verts[i]][A];
+ b0 = f->poly->verts[f->verts[i]][B];
+ a1 = f->poly->verts[f->verts[(i+1) % f->numVerts]][A];
+ b1 = f->poly->verts[f->verts[(i+1) % f->numVerts]][B];
+ da = a1 - a0;
+ db = b1 - b0;
+ a0_2 = a0 * a0; a0_3 = a0_2 * a0; a0_4 = a0_3 * a0;
+ b0_2 = b0 * b0; b0_3 = b0_2 * b0; b0_4 = b0_3 * b0;
+ a1_2 = a1 * a1; a1_3 = a1_2 * a1;
+ b1_2 = b1 * b1; b1_3 = b1_2 * b1;
+
+ C1 = a1 + a0;
+ Ca = a1*C1 + a0_2; Caa = a1*Ca + a0_3; Caaa = a1*Caa + a0_4;
+ Cb = b1*(b1 + b0) + b0_2; Cbb = b1*Cb + b0_3; Cbbb = b1*Cbb + b0_4;
+ Cab = 3*a1_2 + 2*a1*a0 + a0_2; Kab = a1_2 + 2*a1*a0 + 3*a0_2;
+ Caab = a0*Cab + 4*a1_3; Kaab = a1*Kab + 4*a0_3;
+ Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3;
+ Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3;
+
+ P1 += db*C1;
+ Pa += db*Ca;
+ Paa += db*Caa;
+ Paaa += db*Caaa;
+ Pb += da*Cb;
+ Pbb += da*Cbb;
+ Pbbb += da*Cbbb;
+ Pab += db*(b1*Cab + b0*Kab);
+ Paab += db*(b1*Caab + b0*Kaab);
+ Pabb += da*(a1*Cabb + a0*Kabb);
+ }
+
+ P1 /= 2.0;
+ Pa /= 6.0;
+ Paa /= 12.0;
+ Paaa /= 20.0;
+ Pb /= -6.0;
+ Pbb /= -12.0;
+ Pbbb /= -20.0;
+ Pab /= 24.0;
+ Paab /= 60.0;
+ Pabb /= -60.0;
+}
+
+compFaceIntegrals(FACE *f)
+{
+ double *n, w;
+ double k1, k2, k3, k4;
+
+ compProjectionIntegrals(f);
+
+ w = f->w;
+ n = f->norm;
+ k1 = 1 / n[C]; k2 = k1 * k1; k3 = k2 * k1; k4 = k3 * k1;
+
+ Fa = k1 * Pa;
+ Fb = k1 * Pb;
+ Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1);
+
+ Faa = k1 * Paa;
+ Fbb = k1 * Pbb;
+ Fcc = k3 * (SQR(n[A])*Paa + 2*n[A]*n[B]*Pab + SQR(n[B])*Pbb
+ + w*(2*(n[A]*Pa + n[B]*Pb) + w*P1));
+
+ Faaa = k1 * Paaa;
+ Fbbb = k1 * Pbbb;
+ Fccc = -k4 * (CUBE(n[A])*Paaa + 3*SQR(n[A])*n[B]*Paab
+ + 3*n[A]*SQR(n[B])*Pabb + CUBE(n[B])*Pbbb
+ + 3*w*(SQR(n[A])*Paa + 2*n[A]*n[B]*Pab + SQR(n[B])*Pbb)
+ + w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1));
+
+ Faab = k1 * Paab;
+ Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb);
+ Fcca = k3 * (SQR(n[A])*Paaa + 2*n[A]*n[B]*Paab + SQR(n[B])*Pabb
+ + w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa));
+}
+
+void compVolumeIntegrals(POLYHEDRON *p)
+{
+ FACE *f;
+ double nx, ny, nz;
+ int i;
+
+ T0 = T1[X] = T1[Y] = T1[Z]
+ = T2[X] = T2[Y] = T2[Z]
+ = TP[X] = TP[Y] = TP[Z] = 0;
+
+ for (i = 0; i < p->numFaces; i++) {
+
+ f = &p->faces[i];
+
+ nx = fabs(f->norm[X]);
+ ny = fabs(f->norm[Y]);
+ nz = fabs(f->norm[Z]);
+ if (nx > ny && nx > nz) C = X;
+ else C = (ny > nz) ? Y : Z;
+ A = (C + 1) % 3;
+ B = (A + 1) % 3;
+
+ compFaceIntegrals(f);
+
+ T0 += f->norm[X] * ((A == X) ? Fa : ((B == X) ? Fb : Fc));
+
+ T1[A] += f->norm[A] * Faa;
+ T1[B] += f->norm[B] * Fbb;
+ T1[C] += f->norm[C] * Fcc;
+ T2[A] += f->norm[A] * Faaa;
+ T2[B] += f->norm[B] * Fbbb;
+ T2[C] += f->norm[C] * Fccc;
+ TP[A] += f->norm[A] * Faab;
+ TP[B] += f->norm[B] * Fbbc;
+ TP[C] += f->norm[C] * Fcca;
+ }
+
+ T1[X] /= 2; T1[Y] /= 2; T1[Z] /= 2;
+ T2[X] /= 3; T2[Y] /= 3; T2[Z] /= 3;
+ TP[X] /= 2; TP[Y] /= 2; TP[Z] /= 2;
+}
+
+
+/*
+ ============================================================================
+ main
+ ============================================================================
+*/
+
+
+int main(int argc, char *argv[])
+{
+ POLYHEDRON p;
+ double density, mass;
+ double r[3]; /* center of mass */
+ double J[3][3]; /* inertia tensor */
+
+ if (argc != 2) {
+ printf("usage: %s <polyhedron geometry filename>\n", argv[0]);
+/* exit(0);*/
+ return 0;
+ }
+
+ readPolyhedron(argv[1], &p);
+
+ compVolumeIntegrals(&p);
+
+
+ printf("\nT1 = %+20.6f\n\n", T0);
+
+ printf("Tx = %+20.6f\n", T1[X]);
+ printf("Ty = %+20.6f\n", T1[Y]);
+ printf("Tz = %+20.6f\n\n", T1[Z]);
+
+ printf("Txx = %+20.6f\n", T2[X]);
+ printf("Tyy = %+20.6f\n", T2[Y]);
+ printf("Tzz = %+20.6f\n\n", T2[Z]);
+
+ printf("Txy = %+20.6f\n", TP[X]);
+ printf("Tyz = %+20.6f\n", TP[Y]);
+ printf("Tzx = %+20.6f\n\n", TP[Z]);
+
+ density = 1.0; /* assume unit density */
+
+ mass = density * T0;
+
+ /* compute center of mass */
+ r[X] = T1[X] / T0;
+ r[Y] = T1[Y] / T0;
+ r[Z] = T1[Z] / T0;
+
+ /* compute inertia tensor */
+ J[X][X] = density * (T2[Y] + T2[Z]);
+ J[Y][Y] = density * (T2[Z] + T2[X]);
+ J[Z][Z] = density * (T2[X] + T2[Y]);
+ J[X][Y] = J[Y][X] = - density * TP[X];
+ J[Y][Z] = J[Z][Y] = - density * TP[Y];
+ J[Z][X] = J[X][Z] = - density * TP[Z];
+
+ /* translate inertia tensor to center of mass */
+ J[X][X] -= mass * (r[Y]*r[Y] + r[Z]*r[Z]);
+ J[Y][Y] -= mass * (r[Z]*r[Z] + r[X]*r[X]);
+ J[Z][Z] -= mass * (r[X]*r[X] + r[Y]*r[Y]);
+ J[X][Y] = J[Y][X] += mass * r[X] * r[Y];
+ J[Y][Z] = J[Z][Y] += mass * r[Y] * r[Z];
+ J[Z][X] = J[X][Z] += mass * r[Z] * r[X];
+
+ printf("center of mass: (%+12.6f,%+12.6f,%+12.6f)\n\n", r[X], r[Y], r[Z]);
+
+ printf("inertia tensor with origin at c.o.m. :\n");
+ printf("%+15.6f %+15.6f %+15.6f\n", J[X][X], J[X][Y], J[X][Z]);
+ printf("%+15.6f %+15.6f %+15.6f\n", J[Y][X], J[Y][Y], J[Y][Z]);
+ printf("%+15.6f %+15.6f %+15.6f\n\n", J[Z][X], J[Z][Y], J[Z][Z]);
+
+}
=== added directory 'nitsche_method/integration_on_complex_domains/monte_carlo'
=== added directory 'nitsche_method/integration_on_complex_domains/triangulation'
=== added directory 'nitsche_method/overlapping_meshes'
=== added file 'nitsche_method/overlapping_meshes/MeshIntersector.cpp'
--- nitsche_method/overlapping_meshes/MeshIntersector.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/overlapping_meshes/MeshIntersector.cpp 2010-03-16 14:39:09 +0000
@@ -0,0 +1,175 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-01-16 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-01-16
+// back changed: 2010-01-27
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+// =====================================================================================
+
+#include "OverlappingMeshes.h"
+
+#include <dolfin/mesh/BoundaryMesh.h>
+#include <dolfin/mesh/Vertex.h>
+#include <dolfin/mesh/Cell.h>
+
+#include <dolfin/log/dolfin_log.h>
+#include <dolfin/common/NoDeleter.h>
+
+#include <iostream>
+
+using namespace dolfin;
+using dolfin::uint;
+
+
+ MeshIntersection::OverlapData::OverlapData(boost::shared_ptr<const Mesh> mesh)
+: mesh(mesh)
+{
+ intersected_domain = MeshFunction<uint>(*mesh,mesh->topology().dim());
+ intersected_domain = 0;
+}
+
+MeshIntersection::MeshIntersection(const Mesh & mesh_1, const Mesh & mesh_2)
+ :
+ _overlapped_domain(new OverlapData(reference_to_no_delete_pointer(mesh_1)))
+ ,_overlapping_domain(new OverlapData(reference_to_no_delete_pointer(mesh_2)))
+ ,_overlapped_boundary(new OverlapData(boost::shared_ptr<const Mesh>(new BoundaryMesh(mesh_1))))
+ ,_overlapping_boundary(new OverlapData(boost::shared_ptr<const Mesh>(new BoundaryMesh(mesh_2))))
+{}
+
+//MeshIntersection::MeshIntersection(boost::shared_ptr<const Mesh> mesh_1, boost::shared_ptr<const Mesh> mesh_2)
+// :
+// _overlapped_domain(mesh_1)
+// ,_overlapping_domain(mesh_2)
+// _overlapped_boundary(boost::shared_ptr<const Mesh>(new BoundaryMesh(mesh_1)))
+// ,_overlapping_boundary(boost::shared_ptr<const Mesh>(new BoundaryMesh(mesh_2)))
+//{}
+
+
+void MeshIntersection::compute_overlap_map()
+{
+
+ const Mesh & mesh_1 = *(_overlapped_domain->mesh);
+ const Mesh & mesh_2 = *(_overlapping_domain->mesh);
+
+ const Mesh & boundary_1 = *(_overlapped_boundary->mesh);
+ const Mesh & boundary_2 = *(_overlapping_boundary->mesh);
+
+ EntityEntitiesMap & cell_cell_map = _overlapped_domain->entity_entities_map;
+ EntityEntitiesMap & facet_1_cell_2_map = _overlapped_boundary->entity_entities_map;
+ EntityEntitiesMap & facet_2_cell_1_map = _overlapping_boundary->entity_entities_map;
+
+ CellIterator cut_cell(mesh_1);
+ CellIterator cut_facet(boundary_2);
+
+ //Step 1:
+ //Intersect boundary of mesh_2 with mesh_1.
+ //to get the *partially*
+ //intersected cells of mesh_1.
+ //This calculates:
+ // a) *partially* intersected cells of mesh1
+ // b) the exterior facets of mesh_2 which are (part) of the artificial interface.
+ // c) *partially* intersected exterior facets of mesh1 and mesh2.
+
+ for (CellIterator cell(boundary_2); !cell.end(); ++cell)
+ {
+ mesh_1.all_intersected_entities(*cell, facet_2_cell_1_map[cell->index()]);
+
+ if (facet_2_cell_1_map[cell->index()].empty())
+ facet_2_cell_1_map.erase(cell->index());
+ //If not empty add cell index and intersecting cell index to the map.
+ else
+ {
+ //Iterate of intersected cell of mesh1, find the overlapping cells of
+ //mesh 2 and mark cells in mesh1 as partially overlapped.
+ for (EntityListIter cell_iter = facet_2_cell_1_map[cell->index()].begin(); cell_iter != facet_2_cell_1_map[cell->index()].end(); ++cell_iter)
+ {
+ mesh_2.all_intersected_entities(cut_cell[*cell_iter], cell_cell_map[*cell_iter]);
+ _overlapped_domain->intersected_domain[*cell_iter] = 1;
+ }
+
+ //Compute partially overlapped boundary cells of mesh1 and mesh2.
+ //
+ //@remark: Clarify whether it is faster to check first if any and then
+ //if compute indeces or just try to compute immediately and erase if the cell
+ //index container is empty. Rational: We want to avoid a copy operator
+ //(linear). A "any intersection" test should have the same complexity as
+ //a "compute all intersection" if no intersection occurs. If a
+ //intersection occurrs, we want to compute all intersection anyway.
+ //1. Version Compute right away and deleting (delete operation is amortized constant) if empty
+ //2. Version Check first and compute if intersected.
+ //3. Compute and assign if not empty
+ //@remark What is faster? Recompute intersecting cells from mesh2 for the
+ //exterior facets in mesh1 or map facet index to cell index, and assign
+ //their cell set (which might be bigger as the set, which really only
+ //intersects the facet).
+
+ EntityList cut_faces;
+ boundary_1.all_intersected_entities(*cell,cut_faces);
+
+ if (!cut_faces.empty())
+ {
+ _overlapping_boundary->intersected_domain[cell->index()] = 1;
+
+ //Compute for each cut exterior facet in mesh1 the cutting cells in
+ //mesh2, mark facet as partially overlapped.
+ for (EntityListIter face_iter = cut_faces.begin(); face_iter != cut_faces.end(); ++face_iter)
+ {
+ mesh_2.all_intersected_entities(cut_facet[*face_iter],facet_1_cell_2_map[*face_iter]);
+ _overlapped_boundary->intersected_domain[*face_iter] = 1;
+ }
+ }
+ else
+ {
+ _overlapping_boundary->intersected_domain[cell->index()] = 2;
+ }
+ }
+ }
+
+ //Step 2:
+ //Determine all cells of Mesh 1, which are fully overlapped. This is done by
+ //going through the cells, check if they are not partially overlapped and
+ //therefore must then be fully overlapped if any vertex is intersecting
+ //mesh2.
+ for (CellIterator cell(mesh_1); !cell.end(); ++cell)
+ {
+ if (_overlapped_domain->intersected_domain[cell->index()] != 1 && mesh_2.any_intersected_entity(VertexIterator(*cell)->point()) != -1)
+ _overlapped_domain->intersected_domain[cell->index()] = 2;
+ }
+
+ //Step 3:
+ //Determine all cells of the boundary of mesh 1, which are fully overlapped.
+ //Same method as in Step 2.
+ for (CellIterator cell(boundary_1); !cell.end(); ++cell)
+ {
+ if (_overlapped_boundary->intersected_domain[cell->index()] != 1
+ && mesh_2.any_intersected_entity(VertexIterator(*cell)->point()) != -1)
+ _overlapped_boundary->intersected_domain[cell->index()] = 2;
+ }
+}
+
+const MeshFunction<uint> & MeshIntersection::overlapped_domain() const
+{
+ return _overlapped_domain->intersected_domain;
+}
+
+const MeshFunction<uint> & MeshIntersection::overlapping_domain() const
+{
+ return _overlapping_domain->intersected_domain;
+}
+
+const MeshFunction<uint> & MeshIntersection::overlapped_boundary() const
+{
+ return _overlapped_boundary->intersected_domain;
+}
+
+const MeshFunction<uint> & MeshIntersection::overlapping_boundary() const
+{
+ return _overlapping_boundary->intersected_domain;
+}
=== added file 'nitsche_method/overlapping_meshes/MeshIntersector.h'
--- nitsche_method/overlapping_meshes/MeshIntersector.h 1970-01-01 00:00:00 +0000
+++ nitsche_method/overlapping_meshes/MeshIntersector.h 2010-03-16 14:39:09 +0000
@@ -0,0 +1,227 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-01-15 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-01-15
+// Last changed: 2010-03-08
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+//Description: =====================================================================================
+
+#ifndef __MESHINTERSECTOR_H
+#define __MESHINTERSECTOR_H
+
+#include <vector>
+#include <map>
+#include <utility>
+
+#include <boost/shared_ptr.hpp>
+#include <boost/shared_array.hpp>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Gmpz.h>
+#include <CGAL/Homogeneous.h>
+
+#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/IO/Nef_polyhedron_iostream_3.h>
+
+#include <dolfin/common/Array.h>
+#include <dolfin/common/types.h>
+#include <dolfin/mesh/Mesh.h>
+#include <dolfin/mesh/MeshFunction.h>
+#include <dolfin/mesh/Point.h>
+
+#include "PrimitiveTraits.h"
+
+namespace dolfin
+{
+ class Mesh;
+
+ typedef std::vector<uint> EntityList;
+ typedef std::vector<uint>::const_iterator EntityListIter;
+ typedef std::map<uint, EntityList> EntityEntitiesMap;
+ typedef std::map<uint, EntityList>::const_iterator EntityEntitiesMapIter;
+// typedef std::map<const Mesh *,EntityEntitiesMap> MeshCutEntitiesMap;
+
+
+// typedef CGAL::Homogeneous<CGAL::Gmpz> Kernel;
+ typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+ typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3;
+ typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron_3;
+ typedef Nef_polyhedron_3 Polyhedron;
+
+ ///Defines a quadrature rule by a pair of points and weights.
+ typedef std::pair<std::vector<Point>,std::vector<double> > QuadratureRule;
+
+ class MeshIntersection;
+
+
+ ///@todo:
+ ///@todo Need a function which indicates the end of the iteration.
+ ///Think about a better design after the first version is running.
+ template <int dim>
+ class MeshIntersectionIterator
+ {
+ typedef Primitive_Converter<3,CGAL::Exact_predicates_exact_constructions_kernel> PT;
+ typedef PT::Polyhedron Polyhedron;
+
+ public:
+ MeshIntersectionIterator(const MeshIntersection & overlapping_mesh);
+
+ /// Step to the next cut cell and computes the relevant quantities aferwards.
+ MeshIntersectionIterator<dim> & operator++();
+
+ ///Check if we haved reached the end.
+ bool end() { return _map_pos == _map_end_pos; }
+
+ ///Return a quadrature rule, meaning a list of points and weights.
+// Quadrature_Rule quadrature_rule();
+
+ ///Returns a presentation of the cut cell.
+ Polyhedron polyhedron() { return _polyhedron ; }
+
+ ///Returns a polyhedron as the presentation of the union of the overlapping cells.
+ Polyhedron overlapping_polyhedron() { return _overlapping_polyhedron; }
+
+ ///Returns a polyhedron as the presentation of the intersected part.
+ Polyhedron overlapped_polyhedron();
+
+ private:
+
+ void compute_polyhedrons();
+
+ const MeshIntersection & _overlapping_mesh;
+
+ EntityEntitiesMapIter _map_pos;
+ EntityEntitiesMapIter _map_end_pos;
+
+ //Intersected MeshEntity
+ MeshEntityIterator _entity_iter;
+
+ //MeshIterator to iterate over the intersecting cells of the other mesh.
+ MeshEntityIterator _overlapping_entity_iter;
+
+ //Polydreon, which decodes the intersected part of the cell.
+ Polyhedron _polyhedron;
+
+ //Polydreon, which decodes the the polyhedron, which is built up by the
+ //overlapping cells belonging to the other mesh.
+ Polyhedron _overlapping_polyhedron;
+
+ };
+
+ ///This class present a collection of overlapping meshes and provides
+ ///functionality to compute cell-cell, cell-facet overlaps.
+ ///@todo Improve design, it is very unflexible concerning extension to more than 2 meshes.
+ ///e.g. OverlapData should contain a map, which maps the entity_entities_map to the corresponding
+ ///mesh, which the entitylist is refering to.
+ class MeshIntersection {
+
+ ///Helper class to store mesh and corresponding data like intersection maps
+ ///and meshfunctions.
+ struct OverlapData {
+ OverlapData(boost::shared_ptr<const Mesh> mesh);
+ boost::shared_ptr<const Mesh> mesh;
+ MeshFunction<uint> intersected_domain;
+
+ EntityEntitiesMap entity_entities_map;
+
+// MeshCutEntitiesMap entity_entities_map;
+ };
+
+ public:
+
+ ///Constructor takes a list/vector of meshes. The order of meshes defines
+ ///also the "overlapp" priority, i.e. if 2 meshes overlap, the one who
+ //appears later in the list actually covers the later one.
+ MeshIntersection(const Mesh & mesh_1, const Mesh & mesh_2);
+// MeshIntersection(boost::shared_ptr<const Mesh> mesh_1, boost::shared_ptr<const Mesh> mesh_2);
+
+ ///Computes the overlap mapping. Mesh2 overlaps mesh1. Computes (for
+ ///efficient reasons) in addition the boundary overlaps and the artificial
+ ///interface.
+ void compute_overlap_map();
+
+ //Return meshfunctions. Only first test, think about better design, for
+ //example class like OverlapData or suchlike.
+ const MeshFunction<uint> & overlapped_domain() const;
+ const MeshFunction<uint> & overlapping_domain() const;
+ const MeshFunction<uint> & overlapped_boundary() const;
+ const MeshFunction<uint> & overlapping_boundary() const;
+
+ private:
+ template<int dim>
+ friend class MeshIntersectionIterator;
+
+ const boost::shared_ptr<OverlapData> _overlapped_domain;
+ const boost::shared_ptr<OverlapData> _overlapping_domain;
+
+ const boost::shared_ptr<OverlapData> _overlapped_boundary;
+ const boost::shared_ptr<OverlapData> _overlapping_boundary;
+
+ };
+
+template <int dim>
+MeshIntersectionIterator<dim>::MeshIntersectionIterator(const MeshIntersection& overlapping_mesh)
+ : _overlapping_mesh(overlapping_mesh)
+
+{
+ _map_pos = overlapping_mesh._overlapped_domain->entity_entities_map.begin();
+ _map_end_pos = overlapping_mesh._overlapped_domain->entity_entities_map.end();
+
+ ///Oh herregud, please simplify datastructures!!!!
+ _entity_iter = MeshEntityIterator(*(overlapping_mesh._overlapped_domain->mesh),
+ overlapping_mesh._overlapped_domain->mesh->topology().dim());
+ _entity_iter[_map_pos->first];
+
+ _overlapping_entity_iter = MeshEntityIterator(*(overlapping_mesh._overlapping_domain->mesh),
+ overlapping_mesh._overlapping_domain->mesh->topology().dim());
+ _overlapping_entity_iter[_map_pos->second.front()];
+
+}
+
+template <int dim>
+MeshIntersectionIterator<dim> & MeshIntersectionIterator<dim>::operator++()
+{
+ //Update the entity entities mapping.
+ ++_map_pos;
+ //Bit clumsy, abusing operator[] to set
+ _entity_iter[_map_pos->first];
+ _overlapping_entity_iter[_map_pos->second.front()];
+
+ ///@todo Change this. Use lazy initialization.
+ compute_polyhedrons();
+ return *this;
+}
+
+template <int dim>
+void MeshIntersectionIterator<dim>::compute_polyhedrons()
+{
+ Polyhedron cell_polyhedron(PT::datum(_entity_iter[_map_pos->first]));
+ _overlapping_polyhedron.clear();
+ for (EntityListIter index = _map_pos->second.begin(); index != _map_pos->second.end() ; ++index)
+ {
+ _overlapping_polyhedron += PT::datum(_overlapping_entity_iter[*index]);
+ }
+ cell_polyhedron -= _overlapping_polyhedron;
+}
+
+template <int dim>
+typename MeshIntersectionIterator<dim>::Polyhedron MeshIntersectionIterator<dim>::overlapped_polyhedron()
+{
+ Polyhedron cell_polyhedron(PT::datum(_entity_iter[_map_pos->first]));
+ return cell_polyhedron * _overlapping_polyhedron;
+}
+
+
+} //end namespace dolfin
+
+#endif // ----- #ifndef __MESHINTERSECTOR_H -----
=== added file 'nitsche_method/overlapping_meshes/PrimitiveTraits.h'
--- nitsche_method/overlapping_meshes/PrimitiveTraits.h 1970-01-01 00:00:00 +0000
+++ nitsche_method/overlapping_meshes/PrimitiveTraits.h 2010-02-10 03:14:09 +0000
@@ -0,0 +1,69 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-02-02 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-02-02
+// Last changed: 2010-02-03
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+// =====================================================================================
+
+#include <dolfin/mesh/Vertex.h>
+#include <dolfin/mesh/Point.h>
+
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Point_3.h>
+#include <CGAL/Nef_polyhedron_2.h>
+#include <CGAL/Polyhedron_3.h>
+
+namespace dolfin
+{
+
+ template <int dim, typename Kernel > struct Primitive_Converter ;
+
+ template<typename Kernel> struct Primitive_Converter<2,Kernel>
+ {
+ typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
+ typedef CGAL::Nef_polyhedron_2<Kernel> Nef_polyhedron;
+ typedef Nef_polyhedron Datum;
+ typedef typename Kernel::Point_3 Point_3;
+
+ static const int dim = 2;
+
+// static Datum datum(const MeshEntity & entity) {
+
+// return Datum(point);
+// }
+ };
+
+ template<typename Kernel> struct Primitive_Converter<3,Kernel>
+ {
+ typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3;
+ typedef CGAL::Nef_polyhedron_3<Kernel> Polyhedron;
+ typedef Polyhedron Datum;
+ typedef typename Kernel::Point_3 Point_3;
+
+ static const int dim = 3;
+ static Datum datum(const MeshEntity & entity) {
+ VertexIterator v(entity);
+ Point_3 p1(v->point());
+ ++v;
+ Point_3 p2(v->point());
+ ++v;
+ Point_3 p3(v->point());
+ ++v;
+ Point_3 p4(v->point());
+
+ Polyhedron_3 Q;
+ Q.make_tetrahedron(p1,p2,p3,p4);
+
+ return Datum(Q);
+ }
+ };
+
+} //end namespace dolfin.
=== added file 'nitsche_method/overlapping_meshes/SConstruct'
--- nitsche_method/overlapping_meshes/SConstruct 1970-01-01 00:00:00 +0000
+++ nitsche_method/overlapping_meshes/SConstruct 2010-03-16 14:39:09 +0000
@@ -0,0 +1,13 @@
+import os, commands
+
+# Get compiler from pkg-config
+compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
+
+# Create a SCons Environment based on the main os environment
+env = Environment(ENV=os.environ, CXX=compiler)
+
+# Get compiler flags from pkg-config
+env.ParseConfig('pkg-config --cflags --libs dolfin')
+
+#Program
+env.Program('overlapping_meshes',['MeshIntersector.cpp','main.cpp'])
=== added file 'nitsche_method/overlapping_meshes/demo.py'
--- nitsche_method/overlapping_meshes/demo.py 1970-01-01 00:00:00 +0000
+++ nitsche_method/overlapping_meshes/demo.py 2010-02-10 03:14:09 +0000
@@ -0,0 +1,35 @@
+
+__author__ = "Andre Massing (massing@xxxxxxxxx)"
+__date__ = "2010-01-27 "
+__copyright__ = "Copyright (C) 2010 Andre Massing"
+__license__ = "GNU LGPL Version 2.1"
+
+from dolfin import *
+from numpy import *
+
+if not has_cgal():
+ print "DOLFIN must be compiled with CGAL to run this demo."
+ exit(0)
+
+# Create meshes (omega0 overlapped by omega1)
+#mesh_1 = UnitCircle(20)
+mesh_1 = UnitSquare(10, 10)
+
+mesh_2 = UnitSquare(10, 10)
+
+# Access mesh geometry
+x = mesh_2.coordinates()
+
+# Move and scale second mesh
+#x *= 0.5
+x += 0.5
+
+overlap = OverlappingMeshes(mesh_1,mesh_2)
+overlap.compute_overlap_map()
+
+overlapped_domain = overlap.overlapped_domain();
+print overlapped_domain
+
+p = plot(overlapped_domain, rescale=False)
+
+interactive()
=== added file 'nitsche_method/overlapping_meshes/main.cpp'
--- nitsche_method/overlapping_meshes/main.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/overlapping_meshes/main.cpp 2010-03-05 23:36:16 +0000
@@ -0,0 +1,108 @@
+// =====================================================================================
+//
+// Copyright (C) 2010-01-28 André Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// Modified by André Massing, 2010
+//
+// First added: 2010-01-28
+// Last changed: 2010-02-15
+//
+//Author: André Massing (am), massing@xxxxxxxxx
+//Company: Simula Research Laboratory, Fornebu, Norway
+//
+// =====================================================================================
+
+#include <dolfin/mesh/dolfin_mesh.h>
+#include <dolfin/common/types.h>
+#include <dolfin/plot/dolfin_plot.h>
+
+//#include <CGAL/IO/Qt_widget_Nef_3.h>
+//#include <qapplication.h>
+
+//#include <dolfin.h>
+
+#include "OverlappingMeshes.h"
+
+using namespace dolfin;
+using dolfin::uint;
+
+int main ()
+{
+// UnitSquare mesh_1(10,10);
+// UnitCube mesh_1(10,10,10);
+ UnitCube mesh_1(10,10,10);
+
+// UnitSquare mesh_2(10,10,10);
+// UnitCube mesh_2(10,10,10);
+ UnitSphere mesh_2(10);
+
+ MeshGeometry& geometry = mesh_2.geometry();
+// double w = 1;
+ for (VertexIterator v(mesh_2); !v.end(); ++v)
+ {
+ double* x = geometry.x(v->index());
+// x[0] = cos(w)*(x[0]-0.5) - sin(w)*(x[1]-0.5);
+// x[1] = sin(w)*(x[0]-0.5) + cos(w)*(x[1]-0.5);
+ x[0] *= 0.5;
+ x[1] *= 0.5;
+ x[2] *= 0.5;
+
+ x[0] += 0.4;
+ x[1] += 0.4;
+ x[2] += 0.4;
+ }
+
+ MeshIntersection overlap(mesh_1,mesh_2);
+ overlap.compute_overlap_map();
+
+ MeshIntersectionIterator<3> cell(overlap);
+// Nef_polyhedron_3 cut_polyhedron =
+ cell.polyhedron();
+ cell.overlapped_polyhedron();
+ cell.overlapping_polyhedron();
+
+// Nef_polyhedron_3 overlapped_polyhedron = cell.overlapped_polyhedron();
+// Nef_polyhedron_3 overlapping_polyhedron = cell.overlapping_polyhedron();
+
+
+// uint_set cells;
+// mesh_1.all_intersected_entities(BoundaryMesh(mesh_2),cells);
+
+// MeshFunction<uint> intersection(mesh_1, mesh_1.topology().dim());
+// intersection = 0;
+
+// for (uint_set::const_iterator i = cells.begin(); i != cells.end(); i++)
+// intersection[*i] = 1;
+// plot(intersection);
+
+// BoundaryMesh boundary_1(mesh_1);
+// BoundaryMesh boundary_2(mesh_2);
+
+// cells.clear();
+// boundary_1.all_intersected_entities(boundary_2,cells);
+
+// MeshFunction<uint> intersection(boundary_1, boundary_1.topology().dim());
+// intersection = 0;
+
+// for (uint_set::const_iterator i = cells.begin(); i != cells.end(); i++)
+// intersection[*i] = 1;
+// plot(intersection);
+
+ std::cout <<"overlapped_domain:" << std::endl << overlap.overlapped_domain().str(true);
+ plot(overlap.overlapped_domain());
+
+ std::cout <<"overlapped_boundary:" << std::endl << overlap.overlapped_boundary().str(true);
+ plot(overlap.overlapped_boundary());
+
+ std::cout <<"overlapping_boundary :" << std::endl << overlap.overlapping_boundary().str(true);
+ plot(overlap.overlapping_boundary());
+
+ return 0;
+// QApplication a(argc, argv);
+// CGAL::Qt_widget_Nef_3<Nef_polyhedron_3>* w = new CGAL::Qt_widget_Nef_3<Nef_polyhedron_3>(cut_polyhedron);
+// a.setMainWidget(w);
+// w->show();
+// return a.exec();
+
+}
=== added directory 'nitsche_method/polyhedrons'
=== added file 'nitsche_method/polyhedrons/SConstruct'
--- nitsche_method/polyhedrons/SConstruct 1970-01-01 00:00:00 +0000
+++ nitsche_method/polyhedrons/SConstruct 2010-03-25 23:33:08 +0000
@@ -0,0 +1,17 @@
+import os, commands
+
+# Get compiler from pkg-config
+compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
+
+# Create a SCons Environment based on the main os environment
+env = Environment(ENV=os.environ, CXX=compiler)
+env.Append(CPPPATH = '/usr/include/qt4/')
+env.Append(CPPPATH = '/usr/include/qt4/Qt')
+env.Append(CPPPATH = '/usr/include/qt4/QtGui')
+
+env.Append(CXXFLAGS = '-ggdb')
+
+# Get compiler flags from pkg-config
+env.ParseConfig('pkg-config --cflags --libs dolfin')
+
+env.Program(['intersect_polyhedrons.cpp'])
=== added file 'nitsche_method/polyhedrons/intersect_polyhedrons.cpp'
--- nitsche_method/polyhedrons/intersect_polyhedrons.cpp 1970-01-01 00:00:00 +0000
+++ nitsche_method/polyhedrons/intersect_polyhedrons.cpp 2010-02-10 03:14:09 +0000
@@ -0,0 +1,118 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Nef_polyhedron_2.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Gmpz.h>
+#include <CGAL/Homogeneous.h>
+
+#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/IO/Nef_polyhedron_iostream_3.h>
+
+//#include <CGAL/IO/Qt_widget_Nef_3.h>
+//#include <qapplication.h>
+
+#include <iostream>
+
+#include <vector>
+#include <list>
+
+//typedef CGAL::Simple_cartesian<double> K;
+//typedef CGAL::Simple_cartesian<double> K;
+//typedef CGAL::Exact_predicates_inexact_constructions_kernel<CGAL::Gmpz> K;
+
+//typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+typedef CGAL::Polyhedron_3<K> Polyhedron;
+typedef CGAL::Nef_polyhedron_3<K> Nef_polyhedron;
+
+//typedef CGAL::Homogeneous<CGAL::Gmpz> K;
+//typedef CGAL::Nef_polyhedron_3<K> Nef_polyhedron_3;
+
+typedef K::Tetrahedron_3 Tetrahedron;
+typedef K::Triangle_3 Triangle;
+typedef K::Segment_3 Segment;
+typedef K::Point_3 Point;
+
+typedef Polyhedron::Vertex_iterator Vertex_iterator;
+typedef std::list<Triangle>::iterator Tri_iterator;
+
+int main()
+{
+
+ Point a(1.0, 0.0, 0.0);
+ Point b(0.0, 1.0, 0.0);
+ Point c(0.0, 0.0, 1.0);
+ Point d(0.0, 0.0, 0.0);
+
+ Point A(0.1, 0.0, 0.0);
+ Point B(0.0, 0.1, 0.0);
+ Point C(0.0, 0.0, 2.0);
+ Point D(0.0, 0.0, 0.0);
+
+ //Build a point tetrahedron
+ Polyhedron P;
+ P.make_tetrahedron();
+
+ Polyhedron Q;
+ Q.make_tetrahedron(a, b, c, d);
+ Nef_polyhedron NQ(Q);
+
+ Polyhedron R;
+ R.make_tetrahedron(A, B, C, D);
+ Nef_polyhedron NR(R);
+
+ Nef_polyhedron NQ_NR = NQ * NR;
+
+ Polyhedron F;
+ F.make_triangle(A,B,C);
+ Nef_polyhedron NF(F);
+
+ std::cout << "Exact_predicates_exact_constructions_kernel + SNC_indexed_items"
+ << std::endl
+ << " allows efficient handling of input "
+ "using floating point coordinates"
+ << std::endl;
+
+ Nef_polyhedron N = NF * NQ;
+// Nef_polyhedron N = NF;
+ std ::cout <<"Nef_Polyhedron Tetraeder: " <<std::endl;
+ std::cout << NQ ;
+ std ::cout <<"Nef_Polyhedron Triangle: " <<std::endl;
+ std::cout << NF ;
+ std ::cout <<"Nef_Polyhedron Intersection: " <<std::endl;
+ std::cout << N ;
+
+
+// if(N.is_simple()) {
+// N.convert_to_polyhedron(P);
+// std::cout << P;
+// }
+// else {
+// std::cout << N;
+// }
+
+
+// CGAL::set_ascii_mode( std::cout);
+// for (Vertex_iterator v = P.vertices_begin(); v != P.vertices_end(); ++v)
+// std::cout << v->point() << std::endl;
+// for (Vertex_iterator v = Q.vertices_begin(); v != Q.vertices_end(); ++v)
+// std::cout << v->point() << std::endl;
+
+// if (Q.is_closed())
+
+// Nef_polyhedron_3 NF(F.vertices_begin(),F.vertices_end());
+
+// Nef_polyhedron_3 N = NF * NQ;
+// std::cout << N;
+
+ return 0;
+
+// QApplication a(argc, argv);
+// CGAL::Qt_widget_Nef_3<Nef_polyhedron_3>* w = new CGAL::Qt_widget_Nef_3<Nef_polyhedron_3>(NQ_NR);
+// a.setMainWidget(w);
+// w->show();
+// return a.exec();
+}
# Begin bundle
IyBCYXphYXIgcmV2aXNpb24gYnVuZGxlIHY0CiMKQlpoOTFBWSZTWUuYf2UAqZL/lH/2YDJ/////
/////v////8AIAgAYLbePKHuz1mwfd67Wt7tOF9vuyvp8g6pfe5fO5xd727EVSXQfT3b1zts7bN2
599TPbdQ8++bzkPd7qtsU++7u2mvl67stt2HkHtBbffcON7mRRC+hdve5thsACjXvY6Nj5effc23
z29967rcBfb3nuB66Hnn3sGPDF3du7W7TtYbn14Hb77aa+3zxlh7fcvj7ODbq+Iz733xrZ90u+8s
G1dXTzZT1682PHvOdam93mvHXa2X21Wu5vdm7nG2Xbrdw7rEpL6YvsYo+mRfb3DzjWsk9vcYFUnQ
myoR7e47sSVae+I5jqUusDlbrWSHzIAUEHuegda93u9a13YzvY6Oh8Qc+82667vl9n195O973Q++
2q7bVma1dLifD7456r53NBKJV1rB99u+xqfS31LSueuphtV1x3Hpe7jozrzr3L6ds5vffeTqfPcV
e8Zeu5gteTWWq73d6ttmprfeH3ffbtWpdZRZmMzFStarTbsN1bRkCzNL7jna92JXGunCKF0tlCts
rbXXOpYZBbLs7Q3t6A9eV2oYW2bWqqO2u1vfdT3aN732+9jVpRteH1t7XmaG2rVlrbWtk0VQC2H2
7h924jPnnhJECAAII0AEAmk9CnlMymUYTyniIAGhkGRo0AJQQwoCIgE0FR400o9T2kJ6mmyRoep5
QAAAGgAABKYUSaKT0mjNJkaBojTTJNNGj1GJozUD1ANGgGgDJ6gAJNJIQhGTJpoCJsqfpMU9Mmqn
6HqKbT9U0RpvVPU9T1NNGnlHlAekD0gBEkQFGg0maJiaNMpiGjIyBM0mEammnpU/J6mSm2qep6aj
GSaGAiSIIyaEyJppNoTEnqbBEKeU/1U/Snij1Bp6h6j1NqNA0AAHqDTaI9PgoqzssqJGCkiDUFBk
RJFESiCKEgqodgwSooD42AKLeACdoEQAaAZdfSRAWrf/tVb//Lq1mEJCsEREZBYLFRWAlt/o/u/0
aJqMFBZHk1qcODrrLpA1JCMUSISMRkgqQiJIsiMCAkHoAJD76fej972W7VDVtyv9v+0wnKJH8RsV
f4o1IskUwSwREAUAYqoSIiyAiMYiDEUBViiqMYI/20qxT9v5d+oaO9K9e/hbc0Igwn7Z2pJZmGH5
W0LNx5378pP+SXO/2ew7f6/k8650vzaN5j/f+jqIlteUF01xvnL+ysjBf1v/ccZSGnqqotY4rP7N
jSK64PpXCwTCxVKqp3RKmbHQ1r6Nj07bGxqCfTeWtQTCxVlReuORnfqhtlJ66wPKNvzsSfxd6CG9
v6QQxP0cBXuz9TgdiuWrhgp58D9HWikNn/z+js20w6hrbU481quJTq4ydRJcGf+3vxrmlAL5ijn2
2WSN+zyrilcwUKcc4eH9//eFTLW6uwjIhe8TaLjvV/bSqvHSjt5Y3rrv6mmg864CDqJA3+efoMJE
cv7MUV9ub/b/rc0Qrab2q88qgb9+3s/08T589FsKnxtP7Wu7Pl3dnBqOfZ460jGLNkOOfGc1FVPR
S+oYfb3UmhFJUUhPQMhJvrXfqT96YqqqiH30u297wcQrWHuYfx70DQ9HOlIaFj8EmkFmbNhsMgdG
Yq5tg2mc7POztEOXFqHHFy/jmzNkvhVFECqqgrqJwCOb2Z/G+iXcoj3fz3bA+/Ky+Vdc2V82DCg6
JtewTGqW04KLEfn0Ztut40Ool9t9+uurH1Tjb2nidk40yIXH0bCoPL6FoO3W+LENippPlKEYvC2E
WpSqaxwErVh01k3LLnQonLxQ5M9Ts+CQ9ZteHnlz4e3s+51r9reH33NtLnJNLwN2Ot0aZ1XdM9Te
5vfg6Pa+DFrJx716CURZz/ZnGrou1umDnTpon7E705IvG9ns8vd6/e/DZ7ex6d8CkbLLmMmH2wR7
6+49uEH8LsM3hB/oQyzjVPSy4p0NuTU1RS5wYfQOk/xPgX8V/xL2Pm8z1XOzbrnZ8bhR77pUJZU4
p3/D0N60uqYurqcNa6eaw8PrZh8EajK7Zaf3fp+J8VQXp9ahq3t+eGV4r49rEPdgvZj0UKJkReLK
rCBNpNM82HaJt5vKEQcZB7/5tcXAWvUhj66h6g2ZHv3/Gw9PmuTPSX0yrSSGtffdy6b5MmZOMDPZ
f4RhIs1dRv641SMbRGhxM8A4sxAkFYVXtA0ZK4aCZRa7HSHlSUyIgr5U+3cpnlZJPUnLViIbdeGB
kOu4lsiyOjDBA/V1U0cNADhN5FAD0+FUUJ2DqehHgSnCaE6rYI4nR3V0c7DNPTJ5DrLU35XTO7Rx
v5X2bw6melONrBbDoWV7OqzePGo2trkojIdkVsqoOKtglJ6MGjTj/ogIzbbgphLCE5tCwiXAgqcI
kA44tY0AxOu1d364zw1fO7Pb1UwYKaNtlOegxOp5nbEbWu7fhuveHKEwoxZvCJPNPCTsuGZ5eOoY
xaDybGvq6/su9I9aSQBWZbfglALl7XSwg4qmCvby25GoYIdrDcy+kftYb92/Iy1x1DhKyaQrCw+l
zfvz3e44nLfdIph1esq3ZzOTP9m3d7bpabmUKx5PcRJDELnE8X/tLZoQlSMEYQghTFW1RTxOuTw7
EH7w64zzdnc9IfT9O2OWVTdWv4fX+/rSJVNQpXasjNw5lRsfPvf24iP0Vh/Xcf0vzLfTvXth52tt
CrQoHxENKzOuNFkIxgiUYCiiqGzAxmmUSKIk+T6tzNZmB9ds2+UdIVPoXM9zcMgexDgO/KM6bJEv
MoXWrnKXrPxl4yY+NQvocvGam8Tss266ZlO22ewYfwJt42dOfzzr30r167bcEtX1Ocrod/8M2Djs
KLukNIVO8tMiT+F9DJrTSMB6d1MXXxvmHoG0l2eXfo9LDQL52XyaybvWwdtYNV6dVNumtCE1pndC
d7Ces9NDpz5Xg7MwX7GbbU73Z6m9222pqOVFJO5q7HPBXGxAmx8i4dmUGyZDRc6cOOcxKQMcXHGJ
DZvwSgjBU6NS+Vw0ynV+Xd5/z5d0nH3ZgzAnQV6GMA8wVwKYKsFwM9RvH3PCyNNuO4fYnRzsulTh
1E07uzoQfuvqPlf1yb8snH000zj7hu6cqvGrqWlFY03zIpiUuqQRkwt2gxFluC01p016t7ScioBQ
hSmpWmGGDKKdYOYOzv3/h/3u9htBr1A6mhtavRCMjxsv3M6qL6r4OPLjzK+qsq5wekQ2nXHy9lah
dZjGobbz9vU+lOp6MOSTkhuckoEo7NRjDGQrAFYICKSe7bbGr1O9OXLnxv2cdetcHPHDcY0ad7Wd
aq3LNBzYgGV+lU7BFVITrvBhYJIEPrYiOBHjFCprAOMU4wROUXAgbQcYDyg4QHKKbRMo5RArGlHC
JzhcKUR2h+kP1kZPCHgdDHMNbmn2FYT3cEJ3AwEYooiiqsEEkWa/D1dpt12vlep67ee9X0bbPhq0
c5I7b5vsFnId048N+HD5ZciBopVJrRAuCpU3ERqioU00Xq7LG1I5VZfKjqZV1WqbIpWoxpsmWwPF
zmQrUyVSLpXzMO65w9PnWcRSkTVYGaszptLDEM/HGnbY22wuUEePVzwhnQxIerSknawU/DKddTaC
eZvCJhGCScYjUCVKDdKNvD5XURRm7sxJ6MeDeWW2St0dEmRzYIdsmgog3xrxpRQ18FrDUKjGFVJV
l/N+4pUVRt706fMqp7fljb8mxkt5L/hZLQk5PZzsmsS7HDtopUFRFdvVlfC6C/rhqTBga/cPTVNC
+zqoGvmxcY3IJDhgd+yrpuiFWo4Tgq5z11A65cu3BAAF1z0+WgbVVXKkpvEGppGpcVVRjAMHHW6s
kCHw6ZllRMYhCAqjBVUVSRZWNLCDGosrICigsUUkJUAVZRRRISFSUUUQgWDRCKoAVD61ZKhIkBYR
B64iFog2gCyKDdqhlFFYLEGVBEGIFEgKLAFCMGQ7Q5YSix2ysFRINcdCh73b2Zd/frMO+2HA9f17
Xlz5c63W++/PjXHTsiEIB2CqqqhEQiIAKBIqixBWEoohNAkQLMZA5emktAbQVAqDh8mFivq4/i5A
OsM/DNPop33ocAN8XYODUSEg85KwHA0BE+msb7DbFPuSvAful6nBWJAKwhpQQAqBPax6MHfD/PSP
/f2KsSr0O7w7e7v8Ed3hWCvs7T29upk8pnnNAWCbzdqhBqFCfLXBkFDc0YGBBYZqDdoIGISVqsUV
SSBuwN5sWBKRNx31Q5BEMh+ScLhuIjMPR7VXIH/SM/Wzc6yeW9JoJpBZDrQMsGPlmk+W2znWwJsn
1DdTbuKGt2BoQRnaYNYkb01GKZIGQWqOX1IcnfhsaISKvBwWLCYOSbnmwDBhpCig1Bd8JsYFeTNk
kUi5xSSadDNSKtQZBCRGU0nuPqPX/6Pj+s536bECgusxzMnHwF7Y8z6jt/8w+9AdYOwDMMCsMCoH
hngfrM5IcpDpBnIBD2llLxTQmIBoWA/+/22tgbvvd8DpAPWb4vBAkUJAdxe1+vywRh6T/f+YlIYT
bfc70WKYNjB7ev59Dbc5PUhwhSZ7Ify5OL9Y0OobgL/Pchxg8Uc/KprF1LLlZyufgDJYtYzlXouH
7IZ9ldcd4V3Vjlltk1axYVF2A7DDE0h1A8hfdsG69oXzDaSf3n+bwknAenefhzVjsoXFMLEyCrJa
ime7H+W94E1kknE7Q9KeXruw9WZeLRbTc8jHs+9E4ftH5E5s/mv5VqI4cMAHE/UfkSb7QJC5n3sz
SKAdgokQIiiHwJEgPvFA3A8LUOmrFrJJJJsGksPEm/Q3SQNz0Obc98480syUfEPORV6f2pDgGwVv
fTHh36CGScVMEM4hTqajP9RgqhHgtU1Qs5L7F7g8QYKlm2AVncC3oHh3mm5667YnBH7/c93APJq7
urzcb4B2xftIiHqiO+RHRDpB6UzSUlgaSqenpBF5/+DLuUDgsgL9lKlQCRQIQCRWQgqgskFCLBRQ
EZCDICAEhIc4KvMzopWQkWDAPagv/YiDj1OJLv9jr1eHte+SIX/lGQ5AqlJpFUvdGDyhRhoSslRq
0XGrK2YlISbL9AswDx7z6iv7JAItIDhubdcK3fXp4/bij+0a9CVEidpw7bbaXM+xlUSII1kXHHHk
5aNs56LuwtzsY3erFaTBCYjstHQtq2rJqRRdMECpnDC0MkNMVkVaTYYhQp2dWMqnHF4N7aK2HLXO
svTrzOambQtPSNLhZU0lmuLOsd56wK9dS4rd8dVvOlc8PxHsGR5/hmKF61GEjJkXZaQMM9JQlA0k
8aqpRnUpZ11O2++ePgBzSDnrqiQpoVcIWwgC+FbM33u7M4PaFlllKULaEtprzb07jz9uVDwV01vj
TKvXVWlVtLXjq6a2tzU7sFI9rq+l83S7q264DZ2IXuz6WLHFwTW3QZ7S+6cXlFnvnUxT5nTiLFBQ
zSpUW5ecbb5hXVMqY02PDGMsY1w9eWhWgcnFwlp/SjBmXw7eZvWsWvxa7AzyMc8yFJkYzpostqy2
t7vGl2y9mI04+XdIUJ1hwZonikl1zO/GJLMRszODKFLWmHetqioS1wqqGlkamNarRgGe6V6dancS
WrE6dIWcNWzOqQpUs0mWxGVWNi4hc1TYWKYlcrmGEDNtpYVNKguGRjO2cLxWCElQxyxtppRaGRln
e/QbnvH+B5086dyalgsOLQ0OXezsPHj3f2vkPGfdP1NerPX6vq+I9829ff/KJnwcgDEAoPZQeUbA
hcGOAUampZ+cdMhtE3RicJWIwUy01aTBEgYyAYbRNB3aWaHRptljPntAJzYIULDLHlsGEbOzzc3h
9lrFFX60/F3viPUfsf+L9j9lbF56nJXeiHXmuyD6jxrNoxfdR6C7Ll05WK9jT+B1gh25EzmSEGB5
cm/Z/I7BGQtDTzUSdQxVU/933Q1eIfaEjA7y/q04++AZ/VPO9urdNPyPybJkmwBhhEulODSWALTq
ed7JO4OEjy8m3fY6niuaWdnktDc7bbbnd3HHKLDj3oZuJTmBlrLkluBxXFbxYOVmWZyQcORwVqnn
jy6JhxioAwIQZrOQDKWhrqXgSW8SwgoQZ05HbGHsESSZIR4i8ha9BL+JI+LuhfiUXNchu2CSoULs
D+Ac3/dyP3zwPhQCcgmClLBu81s/2M5hPpChvnXw8T534NJMIXnfl16o8EAMnYQon3h9I01QhSnq
1N8LMH/PYKFwoVg5YKU62dfl/PwXH1BNvn7EXCsYppqRpBqEgrh01/MSioxaofqC0JTZfqb7/8ZL
3R5z/jlttuiZkWZLbbb3kyVmS2219QTQBMfXQFsBzkXxMce3Wyvi/DarayRcw1wzDRRRo0c/yHmD
kbPwPZwatOBH/2TbeIUf9kfGDyoJcuO3tsf+AG4BAb/r5EeG26wBdW5oMCzHZxpO8f3hFQtjHbgy
MQXF/BtOnA0cM+oiJEQRIiREiIiAxIxBiREiJERBPj0j2nPasQOGZnwKO85ev09xO8G0FKRgkJRR
QLIlFmYZDju2wmNd7wUB/xjBiDNT0DiSlxvxoQlPO7NbhtkNDPDlLD0CiiiiiwYxjGMRFFFFDmnm
O/z7Hx/f6OrjYs3npmGTJmSSSSQmTG2x48mQmTJkknOeqaA+Tyh54EzMd3qzRzBAuBYqxpS9SVI4
crFQbKiooAlIoSbUiVlVVIjjP6vxu9PrAOGCtpCEHNHaSnXt2HCyTrPo3hnhfDUzq2a/FqWmZ1Zr
JHHPyvWnxuSSSfe8Xbu6BEZ7Pvyeo/V447bAR032AJJ+9CYmZKuEHp++BIrBaMnNlbXA6q14lvFu
h7T7hWC9goo9sC7J7ux9Ga7blMb7tZmV3ZoRYvZRlxeBUM4xIgspqPXwlf4BQUL9WbOA+reQe7un
k0JIhA0ucgPVIkQlYzejbhUYsZmiEvd8WLyfyNtpNtMC5rWaMJk93x6NYW2l1mKVg4KSoJxRmqSm
l+R2kJK+7bVQGGBQI5jrn540aYI8Jw+aQ632zVO6+ijkrk/Kzj08BrH0+zr68NH58/vP/NFJWgtR
RknbCij8+MGlzY7MKQ11S6EJkpfmyvhfadJLax7WeRCFvZW/cP1Byy6I17EjN2DxCSDg4YuKWbeN
7wiWta0LiX3FKQUxPGCebySuLLDNj4CrQPiN4yMYjBpXtDjF0kK6tqkfkipzni9+eby6Oy7s0HrB
1sXl6b3JrZz1xH6r0Cf1WZVQ/gXzftD207oG21cKwo4v7V1fnp5h7CdrU06hjWXVVykyuXcer8sN
2eGDXlbGGtRn1VO0Mi2F1SEo2dWFeMuszjRFrMxwuYla/tbpbNr3zDQJ5Bg33F3DjUqgf9NhzRrR
6Dt4eoDD5MKnJ4edu2iIPPfA2HqAsSwPEQPc6VfBTg4EEF+mzTWVWtk2bJ1xy9pV88t2rjfG0i9k
oJo96g2bIM4BHTj17+T6nvoehEKTgvdVNeGcsgPwptfiDsmfZhY5ZPvIOiZzA4507bViGgVzIvkt
wakIKthBGWT8s2nFngzwZ4M8yIbILGTc8q7C7eNZNyOlpJMirHte8EcUWwbvH8xq7Qr34MXMR534
wnQQ1bZyjHQggyEssXSN1A8iAKLeTMNgYb33YMzcet54X+WXNXi4DF9V72OffV+bf0zUeqpfsuXA
63Jx9h2hk8A/3E2k1+rd2c6EnwGKLERRQgqCERWCkBSCgx2+BDp+r+o7uHi6WQf2nIsyDjPIpfB8
UashZ7yPjDDWUJeZEHr0rwJU+WuUSqgilUxK2z3fCJhZtpd/BGu4HCPwwJLDgsai7zwh6dk1vZzs
Q5Bi7wHSrK2lCpo+6GUth7KVJWgMNsOmMrbCF/mguKJwHiw4wehOy6pH5p7yyo+LblQPsM0uHUXO
JrKSiosWzyiMRrWMNqMRt+Lttttttttttttttttttttttttt22222hbZC2yFtkLbLbbbbbbbW222
2222222222222222222222222222222222222222222222222222222222222222222222222222
2222222222222222222222222222222222222222222227bbbbeO/s1x1ek8eT5Wc06Ont98eYgX
J3IqwZ1Y1WEh6j3TpIiUfZ7qG9I+hQw/i8HWMVwXQ9I5H3fM+i46QXsKrkV6PiayvznQs3PbmiZ5
vIDwzeZuxM/Aa5W0EN2jzrZzcGlppk58AzhYK8iy8/LDePgcS52tu+Wjn4uNhaRwUHU0C1WZ+GvW
BVg+zUzFsr2IZdUFUhVVbK12tkxOuXposySxTS51E/CtODhbTK6MCaebMmDMhnLzGsv/VMuhXdrf
lrz9/O2rZqvhd9BdOwQlWwhLQ3AytXAvuQso7xfaPSJkx2FQqwyBjrhmQTqZaGd6hlE4Rm7WmVGm
IfsiOietv2u3mOeCHtHq7jNYxqhYWiNHlznuO08KI9TWOyHcOd+noDbzaHDcHwNndier6Zj2Lh+r
zNXoWky7w0Bo306m6L3EawxhzXvIwvk1xtmfs97nAqWKZXSGAY23erI0RT2/M73zn1MO3qd+XM8t
E5m23F49PQX1fgq8lVZ0sTDKbuO7w7td8V3AhLcvhqPwchke8K0uIHliEz3GdSMohcGfPHNEZcW5
A3AqeEMoVMzlhakSkgpVLfkYisi6PfmM7oPozaHvhLgq5Y/G/BHt7xOUnYOkAlG1j29MMukA9yp2
9KH6LdVVzWytS8/OFiktKgvIPENkwyYdFtHtj48YEIQgQeEL4+DP18LK4jloxsAMXEr2nC9yCTAx
xC0D4sDDA3AQARInU4G9fIvI2XOF2NWccWJNfqy0q2D40/1RFxHS/GxtxtuTHkRzgD0BzKtBRkyz
aebVnAopUrAD5MUBI5CvGSJzss9FVM0PgVKmWhBTphmMRp2RISExIDEIEVBXQQ2ky8Lcm1bBVNad
qiqDHzetE6hACc4ezlnErw9r+uT19rR51kEy2uFSqd3FN5RWhfxZxorMvTZ42wcZla2vWMq7LbD6
+QoRpaNQVnC5nLrDuZbFYcd0ZkZnFwoxs59DxH0ySj4hW1lbfEbEdEn1s/GNnVGbyX9taXPoZ34y
6OVqMuNg/e/K8flXZuYXr89S29ytqROxu1T42PpS1AtLDSxRa4FTr3jYyuDE+3HA1MdYz4GpYCED
0nJHvAQW95DL6fOMzSQhKwRHl2lGNQkLIGGaGbBq+YeMajckwcc1ZCydl2hoqJBghLh9H9NLJcRz
wGOuV0DedAzi2cYemEjhEMLjWo9gBElBlfs10ezxadhePSJU/xqQ8tIXtF4wTQ5+b2UoweL1ylzr
doBSZuPDpj2Ii32pEX7vA0rtflA/Fnx42tejhS1FkGEh7rZrtuO9k9oSRAmSfGJIKCxTo0ZB/FOl
5At+BMSCbA6YIGZyIe9x16RIEL+mROcEYH0ec/T7ffXtfgm/3pw1x9/6llZ8snt89NvNrc5vWNS+
2y6tlJJl9ackxDRb5NKa+wVK2vWs0ouM6058lpzDq3TfmOx5zzBJ/AgoKz/Aqp4/Kl1bcMzM0B/6
GDZR+9D4e7rzo+07PQTzwvgiEu/xI1o9bvmGhZvg3qF4FwaZOOESnffBWqYUKSTGs5VTGmKjIQ9Z
SNEziRw+9OuRfd+n4GhPk77paWH7PDb6eWxm90yu7iVlVmMzfbfReXvTL7+cz7jpZi2oikWYa1jf
jnnPFFPHxGXvphS5bLiBkIKpEYGpAWAeEk7OYdn8NOw8E15Y3pjn1Objtmrd99tMW96Yqwbg/Nfx
0mtb30nRaWtKROYevRu9NRl3Pd8emPDe3h8/swdy57ddVmtfJci+u+MYr4mTR7bv3dOWPeY3UyvW
5ydakbOdanrnfeTY6cxvnan3F+Zx3OlOu+uY1qOJ13U6rWG1mqip1Tu9TqqiaaqznVNnOMOrPNC1
Is0LqhZoWaFmlqqqqmhZoWaFmhZoqmdnFmhZoXVVTxrVPNZqoqdVU6zhFXBS4KWcqqyZWMOQuiFj
UNCsaWTS0Ygyy4ZZ0sRpolZws0uppqF1STS1hYpZpcZmlmlzlZpafDPUxrGKyuqeaxnOMY0udaKW
iF0QuiF1MrFUsyTTNNLUlYfCvWVqdU1PS1NLU0udLlcZrU1qZqqWaMxJNVitTWZqqxVVVZWsriWh
tEaNRTS2MVhaisvFa1iW0VAtQlVSzNLItTVTVRVRVRVRVay+daqJmqmaiqiqiq4553qTi8W3nnjn
m2219JNa1bvXPHPOrbbdc75ipqpqqqqqqqaqaqFxFZimaCqzjOomn0Q7408Tqq1NVqpHipWI1UzN
Nlyay+dRE51E1nUS+SnxnOtROZV9POcROqqKqqitYp1qqRnxnUTQ3Nbb8cuaKKKKKKLFrmtt+OXP
jm3NZ1E1Kqz4zqJilZ8Z1E1NKKoqiqKoqiqKrJvfH2NJ2z7TzNGDz7BLo9WQwGBTKGAlJL3gAsbd
b3c1NTUnJtTU1K5Vs1ddZXXWV0roaqKqsLVUVWiiqZmrC1VFVoqiiqoiNJkwSYMQuSDRROFog0ZM
MuDA5rNOuiiCiZWiRyijGk0jJpVNJkfCZSMKSlI5CpQqttVVRRVVbaqqvOkFILhMppH08YkahpJe
B6GgmRpMvA8DQNWnVlWih6KKKKqoUooop1oooorRVVVVVVJVVVVVVFUrDUNVgiAicKiIggJyIhv8
C+uGvV/r27Pj07ayzxb8uPQz5rWBATpkQTjEBhhx2wtlIMKyYWTIjAyTBEogkKURlGEwksC4BkMC
mSRYAVAKkKhIREFCBC5LCSskBEINqtYCx+p/mIQfUEaEPVFyYAGhFSQVh/SkPrfa5UYjlCfAa6tx
CBARkB9oJPeRh/GJ4JCnxSk28LIbsDYpSlKH+MA3QD8RkHlj7518bE4UFGpmd1Ws/yDBZo82CVPj
hphhhceUfRBj61UikYb4FFtPmw9IbBlALwvOvB3HECEWVEWzYHx/a9Zfhtz9+XHHohFVVVVySw2Q
9NZXRuywEkmRVVRF6Q5xV2OjzOOOLbSlKUpeAT8c20Zm34+viW9s1opz5sPOqgkNuXS29UOSIqkk
hYkrozWFE1TPq61kxa+zyErwuuK3okz5fvvLq+UvXf/V4Vdl8Hgcxo5EZQ55Vdjm5VNy6nG/+czn
wcTN+JgCqGQtLETv5kOgihpUESH5fbYKAeADB9F8j6H0iILxkXe2aqKFHTfYcdEStKFCE6QWqJky
BM/l/dE5ImOsEfxsnOT09hY0AgZHmSFiczU1MlckIARPP0zrpylKGzq2ttttpbars0wp+kPTA7jV
SFpj0oKD2mbpCCYmIY8wbHHvHqvO0AsP7PD1ZAG0A946W0Im+Ef0MSwwUbybCAvYEmXbau5Af3pz
yz4LvGOJmdp3n5y7du3Ln0H0H6jm3jagnKh/gdNfy6HJg28UA5CwWbznPISjm6zNGoITGwOISbQn
xM0ZrMOLl2LpEvkmno4rSS5PNhc0XRCEPP2ymwGwgIhEZERBgiZ/r0X3EspshH8Srp6B3zAuuq9k
sHrJvRxJimbR7OWvfyycs7mjj6jzmu/T0OZkztdzBnIsRrOAcxyJKZQhBoDsXpxxxMYEAuma1H7F
TTnJXwQFESTgWkjqIAVFdvNiqiBiboP/QRNlIsKj5+XLj51iouu8rgNEsPm3E5mwveDB/Mf3H2HD
x4n0tHljV3RP5Py+VwN1p8WVR+12gwY/jBSbh1D5Tqvtl5f6jS71vLPWKTyj3DbtdOZg0dJRd4Ec
rob9DJebWFvjaOX55rvGfhBV3mjNk+pRPBrW+W3FHLlS1K9vaSdQbA5v/DvAvCqKOnaB0bgzHn6f
Xxp4hr50cXeFjCSJBQVFUL+3LevuPz3+fOlR9fdOf3JrULEPMC1rEkmCEaqVPSuLNkZ/7QXDBjyy
BQO+Hj3Pw4Ss7/ugV6ye0TLD3/J6Pb+HPZeqBv7vQ+Hvqj8i0oKjcTQAVDWPZ0H6S3vKAahDbLDz
ErD1iLERSMMNo4y/G6UxLtu0oAybJP1JuB1k/UGAh5TBXXcQP+HfdcA+UHO3U8VCeAgj2/2yPh8Z
HgyLmSNhibfU+ecUQO7iUE2koDbM4jxcISK5z98ZgT/R+M7eu3rylWiLuwzO84PZj6/KUqZfD39e
c3awaytGD4SCr4vKrrt8IxjGNPj7bq7kTYEQAVdLa74vH5Y0R8/Zl5R2a9rLlir+xy+GELaj81UQ
TfnzhOre7Vrd+W/zaNtsC+XLLHk8NicMs3HOFkPF+kG68L7nrz3vatjwd4cr5w4W6vx0LiwaGd34
2Qvjg2zAwjuLLmC1hBfH0GzAL7/j1juJr4x1oqWkkkooqtXz74m9tqrLbF26fDu1hr1ePjr4yOJQ
8tYxZpytd8B3c7uhuFbM0URCKoqIbgHnIHnBBCezW3p8fV5ajGMePLBcrXJebduKmHE1PPCnbYKf
I5FRpCDNDlO74NdDFlOxqoBxXi10ZC6/X8pHNasmGOZs6hNTbqSQumi8O3M5FdmZSE122Y4YzreE
ZRlAZ+ztoz+GxeMLFsun5ccHwnPAqVVK4G2udTxd3Cz6z1h4+5ylCuBf7Rzs0w7muHPVAo75nDQH
nJnhOuxQzMm/JMwR9Zhni+yfT/k7+b/btn0I/Jv6VKj8E+l5n148++Z+6LAaCwqFEQIMGNylpiEg
w/Eo5elaDt3J3Dv+eByGKdIm9s+4L/mwhD9KaqDP6UNyBqmqRIkgzximySDVQeIHqAnAaRgIwElb
2vLtI8ENQ/h0QsLrDRmsZV+MP7GUFMi1WGAZ5hezT7DmPaYJ+/RaxGGTs1/G6oQiAJzQkUoQ8/zF
NB9B78kP9I3ffchO6LpESnrKCEWEJiVRchB/oDeCS4IGSicvEfJXP0HkLK3QU2+AmRfh1PWnxGhc
xhH794A8QP+gxRqCLIIEAIAlwhrbayZ/vBvhtvD9h+yifnDFPYvtGJF/L+rxPSmC4LU9oPQ06EhI
XHvN/8QenBwDQjGQmfWO6Mj/WAP3aNnn/f5nd3b2uzzfsz+zQ+ECHfptGdfIyn2XrEcbFNBcEWcO
3jgXYooi8g4s7HGZskp1xt/XILM4+X+OisWWYfKFVgKyl/vYuQwwZiHIjVZKtSJR4w1KidnF+XOr
YxwLWYse8Z1VV53Y1BQ/E6Bx/SBcBwAHAsObDMQXRiLCIAXIwlCWAi4AQBoCIfeVKbO58LmGL959
5WOP/3lZtbO9iqMvyyKTNLB7777jy+/83TuV2Y42Ru+b+u4OIU2jIBIGgMhIHxe6Se+EVEuOBu/b
7wmpkSZOGHX7woRKFClFeCJiFgqirBRZvBRYKKAiLD6+fvJDXRfHLLz4m4Al3HH4lnRzdjCCVGEC
QkQ+s3j/ADkQP9iLSrB6FgMAH6h+ogn5mBNQSeAUoTBJS0niHvIeAF8xQQ7uhDDCwb5+ps33QZyD
/QlGMY/5J/H428+6X22iqLOrjkW1dQ37zQ1YUsMipCEOG7cWjJIm3aT0EPSCHpBBCiIlCIFkGIyM
SqGiEOp9CEHq5oPaZJQpbQoeqAoQIDovGTDg7i3UMWmhsNdK5zRMmTlCSxpGxr0HXsPjH1P8cGdm
HSk75290mujkpovVvf+3ttmrqCHMQxHRDMMmybwM0u4gWaSyIGxuGhYiOwMgIGAQIGwDARcy4K5j
YVcQIAJBigbueLtDv7Dlm69uHUpdaoQNzY2XVlK1ROotGJ4k9Uh0EteQhIQIfvgGAbcTpJyPzLci
w2SnrO1Ib8KIoupO7yVTRpWaCHMFjDmRzchIER0vaqq10uHuLu7bQ5cy53kq7hakpozQN41VhGnf
hwP64SG3lH2O3zh9PmAMBOfP5hLIp2qpn86oes2PoGOKvSL0C6pq6r8AiZZk04QyxVHsTnFsHqdQ
5FWA6I8Be8XARDkgPJHUMVLKEQtdILklMXouw6/cH4fw8IGJLj43NRq2lOEAbnCCk4/G2EDi1TPZ
dB2ql3UjDtdFvCk2YbpScW7JR5Ud/WhR4u+V53oEfX28ZWTDZNb+jmFf4Hw4znnhnh3rX27Ex41w
nU+7U4PVU1NqsWgPtCDfBy+qvHZ+a4MGlHRPokHkeeeDIoqqShRVIsBQPSG5EEEEEKVMhShJYRKg
BRERVhGEVJBCCEERKUCUolQ0VFigoAoojIo2SSCT3a2RzxHbRpFUc0ZiKtyZl1C0kIIT4fRhJpIw
UVAQEBAQDhhFkBwq4mGPiODDEZATK5Sr7ATTxHkHmBCgVhANku35uHJjgKQPvBA5xoSBcq5zFMTI
s0IVifenExu5EamL20W1GiPnjWsscplgy0oYFgoCh81k7E97rOHPhi8PA3libOKPgW0La1j++oIs
mSYMs9RoIG311iGTspNE43dVjRtvaw+5h3gyYvAgSPe7uvX7HaSOa9/gjNsf2szMzQTDDvXj/P8R
dfSlLW/fQcNzMwPQ8fMc8iR5kyZQrr+dRGZWep6ka2a0r8EkkgyOG3OG5O1OohDsbAVgIKWggSG5
qB49S7fHerWBcANozXggnPZt3hCDWdcUvO8+rBx3Fre5iT80mSdMMwkyxJrQ6mhhNaY0fpA7VE/Y
oEEGDFFX32CRkeSk6Z7Abv+8/VLE/3/5lxfy/5KQAmRAiMJA5JZyasmuWRFigoMRirLDabs0GwTY
BzNgY4GfQ9aNDQJ6HktyESBKUx0YYxNrHaCCgVorJUUUUSTxPWqaNaU3ZirorCx1Ebx7hVEmNCqG
YbMefiHgeVWQFC6In2CADBFLI8ER+CASAJ/MCJEjFYLRClP/4L2Fu3ogFDxCJhkWYFEHwfWefH57
fsO226D+yvO2CGssLTRk+xN8ec92oGzDhjxf1Rrpl46Tp9W/LfffM33+MJbQtsAtWkltLbAtstoW
0LbbS228ihE69nf6fpPyn7TmL/GczzZ0ttzqSuqdc6t9maOxe3ZNoP5/v1xtCVJ0j/MmpnKtqS97
KvoBqAdh4pv1Ha322XU/WI952xlN4+EqzpwCBYIAq7F+GKPnka5eEWfbDO4YfT0VfJPkMSUNJkbj
/PRRfUBCoLrEeEH9G2+AeZgTZgsAWbashjNss5JCockmIHHHbzNSByZG0NKhjOXKhphNIFUzyrDd
SjeCvKLprwgQNP1W3Wz1rnrhoZ2wvBE/Jng7Y2+vOe53fjFVxtutVecM5ZUqtJEE61pbaGWIpZqq
q7uauajV3VXaCDatiCDRMkoIazU00Xmc7JKfrXfpBtkU/p9oCKaF/0M4mGYQK20YQC2coVkVGGFU
NbdBl1ES5Eys3egXAMwyf2rLCWWAOYNW1UGgCyo5pyETRGAmRgm4ADMCKhAIIAhGlGM/5yZN/Sv/
F4qqqqqqqqqqqqqJIMIa3TMVWMeERVHSpVtHiqaVWM2tVWMcjEcg4qggYuBcVQZwQaqrGaVqqCFU
SqoM46jnISHcBThJ3pVgLgWBRwO+lBEPrgi8Ry/eQNKx6DANO13oZJveDpLUSAYwsIWFVE4AomvR
ERPg/CRjOGbn2cUdBSXP7uSKENpK57sOK51kmG69/AshyxBOQmbcFGmgVcQdr0g6XIvIhlEYGYQQ
z4Og1nY3no0OcVlucKWiqukhDACEtTOsLr9DEtxtYgYFtLzHO4oYksouK2uQmASuhNZGUaJKNRrC
BEcvzT16QNNZzSJ36ZFs0LUuSCZWUwigLrYVuMzlo4isHCnoRuGZDkLNaLo0Vic2iNAKIPz3OwQQ
ukVhuuUs4xv56elwwg3MIiI+eFBtGCQgEQTYUzIlnUGDlqMXXHWmll2ZhtrdEd+iA6KpQnATEMHn
nU2LGPA0LBru/agGgC4wVOCL3HIyIwh1SqhT4F6gQgXUwmikwxNWYF0LAdNIQkmhkM9kn4gKHGCH
8QXgghk5JEiRkSEkIiiigxjGIoKqosYxFHvNpQC7HUGzH3G3EcrSGtltqKMwwSKVELYlSjAomUlG
VRQu1usL5heIqDeVCy4fnyZrnfcLCRhYZs5EPsEXGEbxuZAzBq4ZEBK0SyM41tWG0mECi+IaGho4
ZGbyGZvBkafkKetT2ifgh4KTAyl7Do3U4DRL2bXSBpHvBKGcmqQi9gA4YApNp8y7JJMS3nxdALAF
JvaN41IRGGP1+bHA4FeZPCnt8Hswd4P2AZJ+aZpfsJ96YJdLXk3nBJ17tvXxeDFXltS1NecpvLZt
j+WZ1WGfTDwqzUUNZqXMxM1sfj+R6I6MeUPPnffhEIh2da6WOeWnBm212cqE7AMwSZJRgDiQH/Dr
Wvc/xx7vf79rWqZLdHGuiRqqU1JOpmnu7u1u7urux31FUXdXubwwgMozMMogLiN1qGruLiddQ1YJ
mElbK61qrY22Tt8F4niakDYbmvhOhfH0G+e8n2n3cI9xTacbqqwQ9cH5AymUZtRWRyq24SuFt/Iu
XMrs5JspBkGKQ4Tf4JGOBCE3EtIxsQhK+HsTY6oZ6ur19nZ2eNz3zSGzmIcJKIVJNrxBzgnC8LuL
AhfR0fF8gQIEDvGAYlIpQLht3A18vX0+7H6vQWCwBjFAYKd8UiIxoIIiMk8MiqwRUFEqqlFtFqgi
1FEge+6J17TAYmytRtHOoPZ3MPEwuA78xURWM5tVT0BE+vRD5lM+BNdD5enOXx+W/okiqoqb8jg7
p5hwH0QRioQ6N5UcDE0HEuKi1FDMFWIAJUAWXSDc+313uhgCj89kBsOBhQ8FY9VikPfpodsAzlrU
Lazaw/OC58B2BfZ9AuAebp0JqnPyOxzg18O4LhedxX4+l05QZJItfQ4mozhgYxLCusWCShzBooig
8yBBQFmzxul7D7DonswcccR7L2sFkYVALNU0VWiB1WRjWJcxLAs2PeZTOoOBft6YB2leaXEa4rvz
wYvWXuWKKo7cHkPcjxca+Y8s6lVRu9UIKFhY2CxcaUoF5gkG+wSk6EicIDpO8fgn0EihBTnJRIBC
Fsu1BYTETnOSN1IjGg6I9eouveIfVvv9A+beOsfXLu+vnsfKt8W+3njMXUb521nOzyu+7Vu/fJrh
bWtuHu9+YmvWm26u3oRPYbnv752xJ6wOeYr4zv9OKqKnhwew6hp0PZo6cjJJc9gaKjduRio4gFIE
L1KqEwDsiTFFx0+oiqqQIJMCZKWuZBEIuPMvL+Shwq5DIyPb39U0eBKTtztTxkHzoqqm54naNDme
Jo6BiUYsNN0kjESBTnoDp6Akm4hqbqB1yQCkoxOHB41CeAztz25nVlbPyEvaJd4lubLQwNtqXDd4
cZtGdAbcbnLbAZwjJJJieE2SSg1PD3IJE+hUhYbRWOD74oSIFWQN+tCRBOkDLpDQt0CFtOB4U9va
gPQlsSokjuw41CQjOWxuWjg9hvakGMy3W4c0QkJkSa1uoa80HDpTtod32PUSpbjra6RLVCxj0Zra
1VBmPBRQfKuSkbpAE8OOoL6wOhG7cQ6uj9fOEg2H50c0zkluc9OtQDn6ewyuccQCZiloinNQCUcj
mXArkpFR3WVRzmZ1SlU8xFudf1civvy3DiU4LBNwKoEcWZG5gkrdjOkUJiZU9b5KAECRyOGHUTwB
L5LCHf0lLx+XyoRmTgu3pYrGd2B3B0vrkihEYsRYyREigKAiKCgxBSIioiogKLFkRESIMFIwQRAR
ERJA/AA/GQn8CBm8VixiSST0gdnJ5O46TuO7o251R5aDf6vEWwJr9Zk+pQr6klzkbQbYcLwZ3swd
WwUZPzCUOYKgLODOaz31vhtyOjdJGPMzu+X60CqaJZLnIJKB8o1+qrD1vrr1CBpplEtcaZkIRgIR
qkc/eIQUA9h0b0Xdzi+7EGYNChQP1jCMiEZHB4bZ6+ckuNk16+M3f2fs1hN+LyfFIYj1Uq7OlUqa
Yh2gevr7IqKRYfATxnP+K3p9Gx1tz02mbyJbYkvTt6U75a7aMsU09oiHlOvHj4R6qfTJ4OPO52o7
rkk8XAfR0OJkjmgc6btW7eTUhwhZvwvvtK6wjg94KGtZw7QQka0Q0RLW93cLdXU3dxM5l6uLq/xi
acyRC7VRRZYj4natUt3fG+31HhP1H1FO4e+dh75fUkw+8PVgYv2pvlVMLwO5yxtJfQn153wKJDcB
mpBKIhTvK+FwLePj3Lt7pWqouPwwEGxqX7EUJUrTo6xJNcZN4fp+2xhcGQ8BOIhCQkJCI9lOCvQF
G+9SHy9WWBn2hXW27bdkibDFHux+ITHD4u9MATOIDkIqh4HtMz3/YYnmc+JxueMvm/V99XmYKOpY
tZ1+YENNnK75lbNqJsrZIBhtyHBFMcd4FVj2T9/t9UQKOCjDw2RCcE8QiB6JmMMbAninuNDShZPg
rRQHqZvpVWuFVDx4wmxrSCIOahA+4PAYNICkHQxpFGV8hvxl2aA5SjwVuyadzsjCM2luMlSG+gID
K4L5XYhUIxMC6y3W4ciLFpaG6qrLzsygdohhBej8maNUcdFdkydUK32VNup2Gd05HDrrNMJoVntU
TiAOrYtAlQdAKkhrxGdo6UhDogFAoYhWXLhvtIzA9JxL8ZEmioTB1V27N6nUsZ1gX4s7uyGRUU59
KPAUjaNjsteMO7CMKqncyPyrAM4MGzeT4m5tNINGUUAqxCMyAZDbhatSnHgbh/wSSLdtmrdt7obR
YTYvZLWOgqrphsvb7vXco9kpkI7HnVyHSnK9KqJ6aa+/a9WRMlbhx21O05LfFjZfs9PDrMwM94NC
fvMJIPCBiZMG7/WAkH0SDuSR3siIGB03L4Y68hRJ6C6+LAyp8fj7T4F2RMo1JA/wSiJNjBCtXuDY
Xi5gwoqOt7QgphejAmLEqFbfjUjxlWXkEqBsd7myWjhoesYENnt+PifM5Pieh6FnMm+EXjKC/JAG
aSUwX3JGIj5fFX8fUGF58u4CHPMjybPx7vUs+dAb3+gQe+ObM6J4eXq8IomhIGz5qzAEdohZiYTM
SA+wM3LrseELuyUJwgYSamLY55he1qIYTQDV99NyCAPcEGpju7N+/tdYnCuqiggO3Pv+gb6PpzVy
za39/pzzVO0wXBeKSIJxMzNZuy7lbu7q7uavWo1cXV/cIxq9JtJWxaVitrk0q52WV2yySDtNVmMs
yIGYwBWcVu/vk655polGFC2wMalh1I3vcrnVxPj1Tw5O9lD7jgW4nE3nAF5fXAT5+Pz4wFYdc6Aa
DiSHl07+nCnVKC7O+TBPsY6TUpdvyc1SCKSchuN7S5JKJWRadcfCAW2+HUOZfWF5G8Ks9HVL1fBk
aXxuDXhDainTprXXrpiqhniNtQugC7aVpGjSSYPJjfqnTxnp1A40DoI+M5krwIZPVgOCA8xIrwNz
vNGXUOTgHXjrNtd1w46Zm7UvwtvzkhtWw64YYqueLcQvFWaY9nF7HLJQpb9QloSwWNYNYPhjNkV3
wSLq5e7XHSYVmlVaL9VAnsVemVbGJoQHthnFG0WGsYdkzMSx0iwGVdhNVhMrsmIQVDIA3NIwGDNj
MufCxMX7D7CarCAJgQ/QgcBdogcQdnPXLkdGj1TdRw4W5hxOpoWuiGkBkSQhOaG/kJfSjDcRgpkj
FfikGtirLyFVDIjYiaGNGe0V8lGprnzYyL4hbebmWM8wbFXJsUGMUAo5UqDWzEtNt2WlXxwElnic
CtfOivAXhfncMISd0CmGXXXfepKvcteAHiilEEe4ielQgbDn56NMDEdaZBxCHIptXxGr3N2qqOqM
1erkYX2Kr6xEz593lPjxp9vB1JJnXYewHMeEz3KWbHk3g8uGjxjOP2pB10QjawMDExtCvEpheRET
MaWLQ0IEcw0NexHLeCW2/Yr8wpgPnVHtVp4Dd5vSnNDpr7R6ZZnSuvcV8phC80ljGrwkWBKiafNy
9mQn5iarHyFx8GTZd8x3DMxz3i58RSmFMXq/sE5zxmNi7p9LGjynLsKqiGHx8CX+fe+L2ne9i6VF
JKhVmYSCIuZpre7uFu7qbuzL6iSrur/EIxtkhJggskTbE7CvWrl/pPrOT3nJ7E/D+V1V1+1jDwg0
IRGIzHH4hL58l85hcXFxeLNC4JDI3+fXFNoana5cVmSR8I824ROhSIbFXAt8GWWCIgeoVkaGMoEw
kNKkSMooQg5emg15qMX1EWfLy+6BGNyE32D2UbJYaNJkK8ZN1+gvT8CRe/cTjwbw2TAaTwQayd1s
iG/jsdAEnmNk2MbInBv4iOzN2gTy1HlshKIinZgs+ecpuqV4SqhFAKFzrNnDSorPbEevIV9VgpZe
5Coi/YnbngWrc2NClIFCGojO3JVlxsFGMzYJGqIGtVEDvYggbi8O5Xg2NUpr1Bg1wekGjfxJ5Cih
aY/WXp0xRnUA2hYLO6SG0nlDMqMy+RnUMUBMXmt2tRcuBgBUPoo5hZfpXlogxxfNqUeAizElIqoH
mY9XRpOs/wB32vonobua8nCoeV2p90Dz6v5fPsPBouMV6mW8nY6EsGCDDFjDqOG0ME1xe+CDCthI
kgCWcI5TzxMJX0IF5gFmsDuYLILFkWCwFgsixQUFjEZYYQtcLViJasCWgCR5dqh18yEIDw+IV0Po
IRW5mRyLIGxEx3UircgamdRRI1CRHndG2aDEoYsU4IBnFMSZZPU4Wn9pwH3ImdgiHDcbm3Fkazsy
ZDgg7Sx07JPvEQ/iQBNhDR4OXOjfkrODW6McbsQ0oiMNcuoyc+N7OvifOgWku4C1HzRkl+CQehV4
fJqT6rJ/BNtDt66yG1dPzysZTM7LMwbnEx+R2Ot4pz8v2Erh8oGungLDZsP1uwq85w5qNK+/GYxR
lda+b5HNZ7BkTPwx8zQ0Z+ON75553LvCOjvNkDvo0axjWneLu7vC3d3d3cTnOnu4ur8iMYNQRtTL
sbFCNVZyrYqMn7L+mN5oXH1EGbEAb0YaTxGHlqUiQoATYA+UUytttksRCDYVEF6Y3SbZBPNyBnp3
W7SQqB7seId9hNqJ8ZE98zJFo9MRZaJDAWJLdBPEeI95TRX2d7VjUos4Tzlh7s3JDs8J3C7IFIGa
YVhhgMMgFakF1MLSnEK0kG2pQ2yCLmy0NMWT6Cqq1MzCTEioshtbJa7J13JLBZtY1LUE8WvMnNb9
AiQmIQUpmZhiVKmc6kgKRCUYmBNhNTs48oE04DB24K8G+WxWhkxrTEoYWhEwRUOlDSjOhIgInUbk
lblTIYYU0DpBIREK2MNJqq+kHNMoXaWUvL00xW5gst+OwcFTFg1sKSbnO/Ix58m5066VT0VaQO1v
08a3UHJ4HI4PQO0bSZUHyJt2A9FGehe0aOqM1VUFDULqs9d48Hd3PUMMjE4Y1hrlIj9iQObhrYml
wKbVjpIVhvubQkVFBtaoYwLSgXXL0B42GA6xBwmrArgsAnicAjRVyCJERAGDFD0chITJgYYko8kI
shofaS+zQKmZWRqUFj0NkB4lU70fkH6Fd9AfOj9iAXf7SJsHdf7d+BNvO0aYxzOU6U65UsJjWvD2
jqSkfryKLNrXR79j4viXPxGdxvodvUqgvAU41oSX3xHWohc74nHFuqjrNZbHNlsp4WpGGMkZTJkx
4NlMpSgKoqfMxwbbZxP2c541M78XbNJFasZmlIiZkmZtXvFpd6W7upu4rWIq7i6v8AmHNbOQTJJZ
Qm1xtmm21tmOMQeBfuBU/EqKow4U+VNVQ1QHioXUOB+gA69ZJJJOOGepRtsQYQTTGIXp9R1Q378/
XhmhnRsJ2ibsvdkqF1YQItDcGN/TGbVDWC0U1ZQ5WtdNJRa8fALUJCs9PdoW5XpB1VxHjfgGt4/p
tZnrmmGYM0F+pWmMS2pIno5HO0yywLsgySSxyYLxSvyOTcXjR5XC7nBtyHG5SF+MYg1HInab7oEK
T4dQZWHBb3kq7W1CL4IsZ8p0E4RJqhGhuFsyogOXBMXkHsVTntM8tsWxqp3CaHaFWLEwVczEvlK+
BsZZcLbDJje0HxDPK20mo3hUWb4m0ZOOFVctYi6Arn29VuqcG54GC3PUx53FrKdEKLyIt52OHO0n
hPBg9JyMh1o8kz35K1Qav8winemASMObN7we5RqIPtmRtueJEJhdPKhGc8lmMOZ1MIeKHJVtmPFy
BQvImxheYBXgDXESKnCVpJOONvBxJwcGQzAV8BxVVFUWjhPiiEHhVntOpbg+0E2xHPK/BA2k0Vpt
Zt0fz2Ug8RtscmuSoaDwqpkq8ihxqfUrtEQ6kA2EVpL13brrw1aWRV8xiGfGWkRxkeZaF8wQMNPs
24bgY9bPh1hx6Q05x3oUiCZDtZ0o6P0utlzmHJ3SMZyVGXI4Z1vkHfo2IVSWzn6t87a2nfe7SYah
3uYiDWol3u7stru7u7ubxjM1d1fvEYvJoiSzbYoStXidq2vbVbydH4GVdwztAHNl5lwjFRDAEzVF
Y989aU0CUDIGMwyfIseobMS5ya7pKd7DFLoduAoWXSQCwiNBu7uJv5Z10sZy8rxI5GYSIxHbKzFP
U9YscPmitIICM+wCHUGlkFmZ7GZfad+kPPgQ4XPrszMqC2shI4BwqIBEnrUaIBZ2EWWl++xGMNVI
ssE1ghFEGEdNDIlgPQt1cMIlCYOV5IYV96sUSIUWc7vn0LgS7ySbd6s7853B3CeTSLY1p6Js5qyQ
UyLxY10hpW+uRtOeGJUGMg2MbLCpRCfJ0nAHPr3OyeiRxNwHBmSS+8RKBhThbD1YgxgQLhsmwGJF
FK+0PMSsEgXM8wvNmngTIDFerE2wtkMomhDKsKVVRMuF2qzKxloaznsqEuY6TYR9GKwp0bBadka8
aNIpicCl6Ng+OhjcdOus4PIyeTcM0YF6KYoZ68HfkgUMnad9vFW+lycHtV8qqaKp3KHYb/d6OHhm
XPXcE45nfz7/nh0Mz7N1jKtCqIDqqi5p6j8Jl94va6XKkKLYqzidlzlvjWmfMl5d9OSrqjawDvrh
MqL8vdMYqsbZ333+PF2yMlDk2tstRKyTMystV3V3d3VXecvqKq7q/qEYwajbYkqhdixNYmYra81t
wdxp3khXjktbv3M3nVV0wEyLvg4QmkMtq5zrhxk4pMKyJtquoSUNqoTrtoYpgNOCpppu4nbxLZ0E
kTMQhbFIAwI2lifK7IpUGAqm77M4pRjmr5LD39J48CymzbFodiHjwO/gTKjrty3C8orIggnQ2w4y
MFbaU2reBN14OxqWyEs0Ik2aDDsXgnHuHtlRWTmXl1aXkzIQZYwuE1pviCtydIl6EhMmFcdX44M8
GSEg3LExqcVIzK81Nr7iWs/UWdd4r7ZmEamMRjCjyGGbJ6DR0thMVmD4jMZD26kSBg1Ww2kJMUJD
5GZcWkyAYynvnjWZ0zWFSNpvLIQ9+5G+VLghUUlIUC4J5SzEm8xLiUNTUyQzFxYFssoGxjwK2MpH
JRuO84yS+2Paht5gbyGd9jz2d87mIC5KWOCCZOTYUOMuYOW2bWrGFKzkNCh+VEQN0ODsjflk6bZd
aCDCchz7uAaCNiKwYxrsOzfJgHPMvsbjGR5fc8VgYhTSntQHuR+YmaYdrN1Z63yXeG0I8sarNudV
j21xjSDjlXgiHgAxXObeDiD6Dt6zafDQtztxadukSJ0GThlab1zAQl3lkdpWuA9fW961vjHWd999
/ddtlFwuKzK2O9mpFmdTE0YeipqqqZqritZqIu7v5iKYvbRsUbbbFiaxOtr2vGax2fXwp9LqrOx6
n2feCUgDBIqnBdFRyVlcK3nqB7jr1Q5ct9j0GGojtU9UAmEOu/hWuFz0e0V1dJTQrBGd3K5RDjxX
1rAusaTf0k774F7E1o5HTRgPrHMDyHPezm+xKBM8+hQWvSW6MeW6fxK2ca65DroTXXglkrgc4Xc4
NwxB4S08nh03gc7nkSa8+U8d2bnWVgIfGy8WUz2BiGpo0p5oM5jm2pjO/Ani2SamSJyWbwG2oXjl
9lCCSDxSWUmLcqF4q6whFsr1nAxlhjXj1vwVtoS0csfzUHZW24d6Lw0BB3d6LLZDi14hAXxHZ5Y5
RPoXGejbxjc7V/rEWKBIYS5EeJHAoGxWIsJwhcFGKIoXG9ZwDv2XXRYbmjibw57OHIzr2eh67wXs
/YrHWfCfYChExzLLszMyuDFaGhXaGKoZmRQzi4WpZXm+etNnc57cDm2DXILGwFOqP3uSAXbjodo1
XV8OFuUiD/qh+n14azXaSJwue0jaC7oAsQTEMiZWp06fgzfhDeVLI/0HkfpLQJlD+ZPAEGAkrHSh
WHdcb/uPrPXa+yipyZ2E4jrA555BYqnR+KsOFRURUQsp/zUhc8i9eyXT/9XuHFVVVDY/0cpJhpFk
RixQXdwgfThrQoqqqmi8uDtGHEPTAA7YkJ0nsoAGECGeclFVJY0zBKKAG4Urr3eb0WF5L2E2+jCy
L58qeWArnhSIF61OPRDAgnaLA9AthUaUjJoLxWCwAEXYLqLqFYi4FMckFLi2AUsLiLYUGC5iy2IC
kFNVQc8xYLYED7D0aV2WLzfYw0nrdJ/Bh9phrR5vd816bKu6zki3i8ccLxuq7ReaTmJzGB6mbdCm
umuj0a123+Fzql3ExJ0EDTDmivddnWqmrp04ysML3ng8PnPhH/kxbzOyz0lZAcjLrq+Mg8D4OZFV
ISYuZV9XFOjIAchIV4EO9e2w5b7mfQ16cqHqSux7d3uYkOUxIiLHvhPeuCzY/OwL97Anhb3RP4Y1
lsQOD4c4hzjpC2unX06T31dGDMVb3EPW2CUGL+XLkw9qZlxFDec/Rnw3qI8Xk9GB7nx593tez28r
5xhVi8uMJwie6ZxmPuwjNjEYPwYubNlcxNmZI/3kC/eI1msmt5GpGiQX/Qj3WNZSzv9kYxjH2evs
qrwZ83Xly4ududpgDPR8eOPdbV393L1uQk5jmN4dYOEKY+Tq322WfLvhy1vYJl/eZuheXt2PbVpI
oPP0Z7g5XEFCgCbCEhNre4XlmLz2cl8IfnJh1dY9dlcPm2EDzgB34PxY1NzWNUOCJ96/nddXDnOi
e7W+rmwgi2AyNcH9jdULLorg0oP0dDjQuwZoWJpB559sBcj+jsqsVimzBAYSW/mi8vbLbh2ecSJE
jHz8/OtBckhVdD373R9szNjPk5Njix3jFW5rwnx9/V+Hxe8bHt68DhDtxstn7HdCDaoc8O9vBRLY
EC9jcYZGje/bjza5c+b82va2bqfI9IM2jKbCZ3TjKSRfPDsjaEh4eBTheuUONvXop8GCufsyv29r
Wbz44RwYatVdsSQuTnn9vMrCkgc7e/v0sZ0MUTOMOmcYcx7RgleO5Rju64zjhk4jSU6W3ENTDQzS
8xhciGskN/HI/N+b56ZA/EYwI7WSrHs8yjETY7EdoUjnfRmFg3qmwcGBs14zK2Bg/BgfkArkQYxY
REX2+78/nw/yw6g9WpAxAjEhJwqsCsiLEikio6QEN7EUM4gO2UJ0ekAMd7QQDHE5/PQB1dl5iKqq
9QQ4JDJ3D54H9cILggDhtiYFhaH0wQ2wMhTAoo3gFC1o5G+ag2gS+ODRiiLs0sq/ihBD/MoT+oMh
EP6pVhDQzmm0EmWpBVESJGKRJJLtSQqB+9Fh/WEYuM/07tKCWIJDEQKUoEgMJQkVKISIkiLds0EQ
irIKBeUQFfuRPeRIxjGMYxjGMYxjGMflrxkVLIkRIiRjGMYxjGMYxjGPZmODeEz1/OZk1pntKRk8
eOIgrkGRSMiaGXHcoyqRkTQy47lGW4O1lbSUalHJhTHMpmpcTBolY4WUsqM1NFNOgVFVGGtU06SG
JghRiMuQZUXFay0kSIQdyjLcgpYy0mS4xl6OIhoesoy0drLXEQg7lGWu1lTQyu5RlUIyV0rKXISt
lUEUBMrx4xloO1lUIxJjGQo6V4xjKUdK8YxlLIllGyGKEGUUIyqEGUoOlYiDhBlG4x4xjIXIpSDK
KEGQeNIYoQZCjpWTG4MZRQjKoQZSulZMY4MZShJRsQmMZR2UjJjHMYyqEGQrpWRQgyldKyNDhBlL
klGyRwwpo1MxSxqYNErrKYxalNS4iDGkMg4iDGkMgOIhHEQg6VlK6UZS5hAiRBMQyGDiIJoY2kMh
g4iEMHEQg4iCaGTAcRCGEmQhAlxjKpkGRjjZE0MuDuUZVIyJoYpkGRlyWMtcRBjuMplNAQ/1WH1M
MKUjIIWwpGSWIsES2lsspSyWCUthApzgAlJFERFVghGAlhCpFgCJIosBYoKCqKlQRisQkRRkBGMR
SKpIpGIMQRVUWCH/XLcwQhaACJsAhg/IP/YP1eebLoyDaLovA8H/Pi9gvOvfEcoOQ1nVWL4W0ta6
ohICLl7asQEWq7jMV34pqQqOIimrUk+n35Mg7BcgQgYcWnGiOeWQGHua6nALFN3PQGy4lnDCGAf6
0WN4cKCUE8f5CwVqIqM26g8UePy+QENnIdiqoRcf3vvHmobfwhInw2a3iT1/WB+0P939z/sFJERT
8M4AXfnf0T3N7gEf+A/3gDk/n4A+irTjAKeTaFxbIh/EfxM7EB3bpoN1/i9grvkIAMUkAGQCH0xq
GMBaAN8x5Ho3If7bQB3Y7DeCgicyf2CBRk9l7km9p6Usjk9F1bWFkYbZPrTW2Fe6yXTUxtxD/VhQ
eVAYN40kFd4hE0Fgn/AghoTUv+yixL37y1rn6s2ziXHo2uFlhyoH7tqZ4a8lrJcaq96G6TEVjMbt
kCtnjmEmtSYUxu14UEkk/n9InDA6Bof5IAfhkNgFOeZgShZRZm7AzED1nYDSQqGrikBUHGAj3Gbh
BuYaUapZ8Ce0AuOvb8dBY5mHcfiHm/zqz4P3Q7giWL/jDvu5Jbng+Ae+yfpdHfGjck1zwHACrJq3
8+G8zm4XrDpInS2fLOnlKuF6eRiHPpB3KwzHRErpvY8YDhuMaKpfkWKIr7nsTL7ZEpRkQG51HOB7
U910wHPOs2OmBhhfAgN2XvQ3Y4YFXwcCyPNgP19tcgIjmTqgdPBBUt47OFgfgOVH+vlNMANy8CVs
2PIsYseh9HL/jlQt15oeRO0ST2sTroHkleY0NKZSSlgk2AEVCGnh9hr/5c4favNQtEZIDH+z1+xV
gqLCpsoeKW0GB9gHkLngnoPuhMTsQvEG+neeMJjvZ+Y3xxAmn5FyGBhIkA3kE9rPXCcNbPkq6QL1
QvoMUCjqOlGAInKej4XyAHpPjT4HRJ7LJzCGQhPL7JWAsqqr4eE9tL6rYImCjMZmOEckWZd54/2c
79k1+5V/H5And8p8XnqZn6lIAvsgvNekoQ+ED0ELwG0Q1Ur5KBKyRPD4PmW/u5hlw4WGd/kO48Tu
RRlh6ygF8Zrs2gajGAiAjAREiIwFWRViJBVRAVYCIkRGAqyKokEE2sPb+VOTOA+z7YTJ8udBxhRG
bqk6kDR0ehGgu4u+YWyHkaA/ygbqRh5/5eu2/r1NmeJ6aB9byiKqqqiIooqqIwGagU5b0QWMk17L
zFHG41MsDTOSRA5ScgzqBjKyFklTgJ+X2+p9y9LBggTqfJxTq7uBxKNZwWNznQU84XfKfRjSGYt3
PoQvLek9od8jAsOClWvhfG3r8duhrWH2w9r0RKZN2bt8tjWpS2UtlLZS2UtlLSWMrLNWZbMjMZRK
lLZS2UtlLZDyJ6YTvNtBqe2fZJHMCmAT3AG4FY6btm67zyKkdjCMiaaGlmmcyKkiiSaYGGYGEZoa
Vl5zIqRIyMkMn6gdHeGpmqs0KkIUMIyGmhpZpnMipEisRN6HhDwV8Zrm7RRGA/49oYCBOLCwlkqL
wqIlsfIN5Ts5LYQ5wNFRRhCIooo4s1qQ30Bv9smHP7fvR3OTD1cyqU6izts/4TACxiCTnIyqqqqq
iiiiqqqidm/UBhOUhQ6uodG9pvNg/2LnKJ8hc+gR8QZwjY4IliHmF7XuS5whxB9ycJEHPkmJ9ou/
VpO3sQO9to7vQ8D2zQ4lNqea9yjWMZ5jgnRyGuOra0txPb9GsRchj7GGb3p28Yl9BsBxI2u0BaWL
NU9qBLPsqQogdEyTWPN3blotygHYmeBj3+F2rUmHqSYl7g7RwKgYVCDxBlS9spiMeGgZkMy5oljo
LIBzq0VGQHtT2JgWQBXpuUm3OGQ2lk7idpNk5UX0DBjMC0rL25hBmkc1BeexVJBohdgGYdEYLcFl
yVXfeCyPR9Q3veTds+Q8JkoGm8hfteOTFoMjIuEYkB/0CKAOyMCXCjyOebhLrJEBKxg+gupcBfYg
B1DgnJOyn6lwRG9+ShsolQSiwav7js9H8lvw0YJ8/truYf1S/E4F2wzF/2YcbhuMuA0ocC4OBJlk
ShrjQwNNqmxXaSTG0mazPARKk/f/g+CH8XuLr/3QDtpIvTp0h4c3iyIMjOnUq9Rbz0Jvz588OZuq
pEVRUVSUkKGrIpqppzKFWlOkqLnS4CWBZw9xWW63i3ra3DYLeMk7YE4m1A4YLdlxWW6NXWnDQJEe
8DWt1rYVgt2WIa1ar3wtha9ijCmsUJtF3jyDRE4QXzEMKn9BsQkJC/JofiM9dH0t4RT7TnwdMNs9
35/7j8WimeCPojyRClKNCHL93X/w/n/1w/4z/ppZZY7vhhhOyuSkpKlKpGxsbbbW0tttRo20rSrb
S3kcjkFqszbbNbZ5/I7SzebzeWWWcjkci1bau2Zjg8iMh4s7u7u7u+KVnr/hSboYgYoGwnf/yxwC
AxgzYYRxrJSlkSKFBhk6v9aqqs9Ntt9EGqBmQQDgcDY5Bre/BKqSTU5WtVzA5tODxvbTQkFU9GS9
m3vyzPAr7HBD8lHZwdnR4Oyjfgrjjwzu+trUp4s0/ThpG3sfVSXEhuDrcpENwunM1NfBevt89q20
tL+4P4ifp1+X+r/bbX+kikYkikkJ0PgOjpGKQwwKSyRIlJaMFItMWMzbUH6opU2hsqAXRNyywRkY
wcJoYscMSIQh98sCxb8JkCBZmZCyyhEAoBsG1n3gwNkD5CSSb5unGF1qLeMC7LSXRDWSBuDLwMCj
AX2JxtwpPb4U6D5ALW+WpUlZP6jUaVC4JWFi4jAzKRAxq96D9AQAiJc5Y3bgJE9BNyEh7QJ/y8lR
FVUViIq4TlzaKgE6h5jnOc6AVdcRZvK0L06HrbUMhcmGBKLAGDG5SgwxWbpgwswpqFlDho/r33Rb
5/yLVOJQOOQiSSBIiqqknQ2AJ6Z+Z+bFkUVjJYhBkBOIp/qJycx5AQx3yGbRR/QdE4Ox+80VtcS9
LoLXQpYpULlIId3Hs5VbhIiggxYMGIQGMghDckaIjGRgsSEEgwihEgkid32zeABoHXQnM4Ce1pYg
Uoh82E6gOV33LRpC1kmgDrBYQfxQrRQuqD0xCQkWd9k16sOrLnSFRLXsGafIkIT+vQDowB1PqOB4
fAUm+7Hx6KGafb2bD9/Z90t/OfqIf43wvUvkQwL1YmsbKw5MOZIdEPMGBarGKqqqqoqI3+z9+HbH
3qgxBBvpd+NtfsCUosOuJgkggXwQcaojEVDILisFXJstQQYwIpA9tSCjkgfQeMkDIJgAaFq8cQHF
kQBsMbtFIA7gkih9Yds3TQIKs5SAJDJAEChTbqEKULr2phmui5xpzcqIOdwoRT9hycSWgmBtaxdE
Hj5ekKLDDuhU7/zaby9wHsQx/KhMDE/RoohQwKdyXSw3YHRqG4GC5k34boFokgsjCwUJRwcMnu6T
/T3aHq8oarqA2mUDyfOLOztGOoLVEUpphTTGc0d1S0DQH0OiqoQxBiecTxK9H8PX8bp7Er9e3h+H
n+cxzAdZrrazwSI1EqJUHzJqfU6/UFeZh3m2au4mutHZFOk2EjTvNuHjtW08cOZkegUHLoEKG5pe
KaAEwsEZ56RZkRxv2EIiaH9WbG4XmJ/g2JMVyDaauuob7kN6ozxMR5EjA4lV40MSlwLyIcuMHftc
Q5OTiAy0IxsEHG/gNgggU5N99R0cA/J3vsOc5BRg6IwDhviTwUHQavpa8EmcEn8aIiXO4prScMWx
dkpwVxD35YTj9h/XjFm2+9LsYDHJuX4nUIr8Gw5gs4lRSjpl14z0Hp8iSXpXqel1YdY3HES7AlMW
mY1B3KVmQuRcXAwclMnad4dDvPVynLi8sXMcYXEbbZWP92MaZ5QmX0dPN49Qh9nV1dV5g2FKkFLC
2JkH15kRxjDQ5zO82HihY8Y7YHORsbXX32tCEOYRIsmaI48iBK28s9APvBcy4E3sBBVpadLRzQtK
iBWebIZCg+EZVAUJAQhEIxQHoPvV3YQBIX7fevkIF7y0qOo6ksM8cu7x6VdmVhaXane69EJN7EWk
AoiQgLBf8vUy1NSjCJPYKLX9ck3p2bja8dySkOT1FX4TvK+FpwvJ2w2iUKFRYbmobt8Tbr5Ea5Rn
pAPAqDc9IFCFoCuqYbgPPAUP/ggFQ+LyOQxQ4scBRxO5QvWyfl5HMZhTB10ucCbHsEdbAHG3X+0v
AKHi2LBYMH/SfeMQSKQUikWZOBqUP5CP6QNT3jELGr2ZLnf+qQcV5RH4d1zWPGgAeMsQxEPMDNoI
bEd9bdH7KLMOmkQDya94vvTH4Hfh2pco4JCFmwtlHwic4A7C5CIXA4KDmllAzUQfP5+fXli8AUt7
jmnhOPdDQHzErV48miVBiRhCLVJIEoAjSFJEdcHUv2QADKBiQQX/eX5N7knw3PtwSs6BrlZDkT85
KLinbsAi9Lhf8DwJAgC2b/NzAcMuZhxKmrbtV4fYG8ZSRmLFxMh2na4080BfuQCAl1VaYMCKmPOS
EYSB1lDwy4TFK3qcVOcgUQjCBsMtD4TEEcAVQsdxl7EWHH3m3f4c9IFMCFK7Q/oQkiwIQgRrf0NR
4zQO+1WhO0KTPpoyedycIhllTmhjCpShqaHWC8RMz8J9CzkWHaxUd8QmVkJRsWywyMRC6AfzTI1K
UaghAYCGLAKijhxvCkbCi96oc/nB4/QEfq4A3RGGgecHvlBDQ8ePivsZEkHelFHpaIkGhg0koexo
osUWGOF0q7KCReRMUssJpwuPzfYfm+z8Cr7iXnXcfILUQwRVBr2E0SLuPBzEUoohwJRWDQUEMGNE
x+wNEU8nvjdPtPoV8f00L9Du0NTuOh3H0y23KHA4Vw4WspI3gFQeY0E0NQunCkTEmmlLcYOhHjeI
KuNz2nbRK657oWkAZj8+hcx0EI5wkYjn8JjpquBoXBnMoG8DiVLe3JhkLAwKEpY8UwxdZMyCiuxt
LT7vpZyvCmWUCwwjiZZUsENUjZ9dipzPO5IujLZ9jbCbXmY+lDAzLsQsoTLwutC0oxmhwgXk/ejA
3NzG7EyyV28R8oTCFAmq65q+MnMXy7XwcOiG23JnaiShr8Nuc7lJlPOxyFF7p4SRshOTk6az4wPu
Z3nre4kktnttqLGs5rmR2WZFu5llIG42kZnDjlx0XA4bVmfuASFxYXzCCOvwcfIeOEhCmEZJCYiK
Uh7R4xiRKTbynU+pT1uOQrBvrMmOLqGXUvftiprYIKUCKbIY7Bw7LTicudfG5Paz3OPdF2gQiNyq
JInIpCEGduQQTVBV1OXq+9zmdc7+uIZFhZMxT1MFC40DrUoNpdaUsVWOO2JZkWfXxAZgGTIZCbDm
c6N8sZ64ji57A4cTU6B6AiHcsSLA2jmL2pIkN9IiNi5jbRYYuccg5YZoSLGEbqhz7DzIpxKPf+g4
P1v0Ip0Khs5aFQlBSmFNAfSq3QHe7Z5Y/4EDvfrgUEDQHxsEBYWL0FPvwO0X1C48DUMAfEXcNSw8
TSD4ROYL9ns8a1+0EXkgIyAg+Z6ieKFaaGa2Smg3qB+EvqHD2GXlO3AwhNE+VfhyXffXgOKTWSqn
Apsbc3LeGS+HASgNRw7QQeta009T3LYJEJHvKzBQmdnAHDmRe2Ndp9UCZtGBzn7jOv5WvGQ6OjfC
a657FJDbRRQx3RsqYGNxcWRel95ILxfgsDAvMTiVluuUTwzNqy8yoeOecbeN1PPDFI3ext66JUO9
hzBHcX3EG4T6Fwh5PVsjTeMnJAOcbIqKUcEYpsEHoa5DY3LksChnWXUGxxGl1CsFtmFkPr64nITR
rYWQy0LnDA7d5zM+Ce7b3H8gQ/HYbJPcsMcSANWboue+Bv16EOvZtTYq1UTrUnVRRJBxKmRiye86
01xUVqRNMyKhjgNd/RToObl1kmWJYFCrn0DORg/jj8RjIqQKaoBgHPXrXQgZWOnC7kVIFWISHEgd
CHAAV4UxC+819zjnZhqYu05kPCxh+1njJ6aKO/8hMQ/1blJEIGWLYsIfHR7xYD5J3Iwue1aD3pO+
1o3IkFn/XhkkMISHt9/IRAX4JWMSSMRFSEFCBCQjFPHSqk3dvYbiukxudR1eTtxMA+oOfLt6i8j6
wYnpm0ZFkCV8rB2iw3AB2R87E0IZ23KFK5umg6MCOA4zfF3s7MikAUcBaH6/noJGJ4OUU+IRE5+G
nBVT9JAM/mq7UvAvcCx1IEh96wiv0oH3kmodu2g+kQMD+0KAg9kIiO7DyssnWyc/Q20pSFq0djYS
wK7QTuAE6/GyHrYgw+UAp2AeeFNCWbpmBZgMswRhVSlCydkkkOyLBARg0ZolKFB7D91BwxdJwEtD
GFiHhAns9BhYk0BXWVRSxitEaIBXcbgE6RAIkUc6+F6WYvE9DsPk0r/BV8adExIwSC1ESlCgYoQQ
J5QKEdowIMCAhIJECKEIAsbqF7KVAFJEocNsC7iNuQHN/i/y28BLtDmSCz5PRlBmYIpEiG/EK8/v
OoBt2zvnx2q1eZsqVznrqmSFiqcSYZEcm7gqhoskPsQw2GMG6skv+52SZqnrRcfC14TqSGAKcDKi
6sShWLmBcMLYTi1ByhKMOW0gXCMhu+Xy2+G757xNxu9JUiJJIgdEgG1d4YeLBJgwZDsTZwcXzUD5
7hjgs+NCDkMnCvQMTElYfCxH0D1EnI232VbwvScvfx92h7bDMzKBIayd4aqWs0V9DFgaBME1paPk
Z5hqIsLW8cTFDoK7JDoXZ6MFY7Ojo7eU7uDc2Dg+QiKPe42BMZRnffKSZWGJDA2WTEcCSopJ7LRr
syMMyMVo4QY4oEh8J22GZE1stiWlCtFxxEgXZ0ObOdm4oTQhkKRSzboRK2wHYfmRPgl6JxG3XTiu
aKhyZ8fEZamaPXBXsEFOKQE5xRWwPQiBwFYKlxaAAoRpiIdIBkfDn4v+MAWazJmvQR1h1LQuJaSJ
VGxaudxUFiM80IUxdntYf5wZn/j9dFA4cHbh2kBIv44HcLM3W52B2E/wEV1HcWqwpWc+vn1106hK
nX18b9AU6NyrPIp5/QooYdMHPgwkJJKAAdkUOCgbSxlJeiAF4ic73ZukVwL9H0O/gDk1XNZdxucu
ZTUtI1BGALiGxJTEgGTIRsv6QCsZgt19+xE9Ux0CbDgvBtQ95eC54RuJotVjlbZiqDlzOddTmVYv
14BR+lVMP4RAYQHiBoVoRQioHGPtAUpG570iH4heiOPomew12XNfRqOqhHlFZiYmf7QA8ApLAKHn
UNcBwQHRA+zv3jrELCym0UPMWQ1onZAMYAP2ir/GlLa222lpaq21bL7aFyltgVa2aAJeROsnWSEn
mAPuM6pAJAL8ZAInwcNxKRmmVPikmuYF0NQXqgNQX5+R19Xw/DIfA3p1uezfzz+vnv0DOS8KmMdO
VDECMxr2GT84HqJb/UrTCCbCvgtCqDFBDxGAFzy9mv76q8mRW5N2xRvI2BjIYpcdACr9mBzvS4rA
gzIOBSMyz9G0aMTb0npuV+uz88cI4ooyLn8jHHOtrb9eilFXm6ZuvjXlxFPEjSPEVrPC6jkMVRm0
HVQVUAEsYt0xELRplPGkRo7aefDSU0gIgNpzD4w+2yCg5pdqMJOlNmXVIqWqq8MBUFit562ern0V
OLoMh18uRo3HPT14YMFI9VOQWbXKSFJMiM0DmCfSh3Bst2BrYKpNMaXKVMjrJlVFF0GF4vjG+QiG
NKb9ZwQo2RRN12t9LSImtDKooLCif7esZpaajb6FrQ+R8/AULpz48h3lzyavRTYi3lRI3bDTgXCh
LVVtaQuvEihtM4I1yWkWJVmL5jj2DRQeKG+vfSq+Hwdpb8J9Qw/RJdG2HQKwrBJULugzDZLrwJw1
sdrnpJe6AKIKEmQ1FFBZPNYDYssYAWYgJwvr2O1QeE8Qa8+Hsw4P+D07fecjlQ/IPunXP21TsaKl
aqqxjUokqiRhKhCSM7DuPlt6uf5vf1+Pj6Pv569/y3T8WvFxz4eDbiaHlMw1KPp2vke5aJGSSSSM
IREFFFVRRfHpyOvoeXQ6zrPb5ezYPKcTf64dRmxDund3RncFSCCRBUVRRVWZmZmZiJZl0Ncolgab
rQwkDBp2D2EDY4EYqnZ7hfch4sVfbds9plRmauvXdI+3fJsjba9IRg7s7KqXSqcl5MKFUn1jSt9j
MxOL7aA+6SqOO0ypEyrfiGImXGxUSLAoOgE2SSVtXF7AXc8FDdiHJXJya5RDHUn57yYDxPR2NsGe
Ak7fwOdHbh1SaPA8Jlk56ncx3WNpOyshZ0az8nO6nacsC0McG4WIczJ3WsJYhVGLtaHzwSIxBixk
CMQYxiMYMYxgJGMYkFABQIjHQBJ81Q3gsgp6AA9CAHJiRRVUVCVQjOYVhyyduQO/TYTvisbwbDcy
04RchFKeJ59xUV1tIsgixB1KseYwYuFBFzRhq1VQpCw4Kg2TzmiXUSFlIFCrQkBJBCgDHNAbzxUD
Aigkk4hIzY5dYOcmRxCPLiT8vVi5GP7BxT2HSW3wWgceILnCc2AJ6aqxc9odvcYTPt8Ch3h4eJyN
gVy/SdYPZh06GF1+ihjfl3FV6eAMr/NIOYuvoAfnHAt1r9gxeogyIpQEZ4qBbJ5CpBLkUIqEVtE7
kIkAf+IjQUNkgcipO3ubxZY4OqF0IEwXYzQpV4zXsGyrHq7rbPl9u0xPGA2cwKinag28M79MjwCK
hnz9ucp4OoWKwfBV1S6juPQ6YB5nT0mcDuITQsRJWguWbgHMiJJAS3FCwwS6JyC+94dpkoppbtw9
7aoMIBy7yhps5xVAfiLZmHPGbg2G4xKTu8qXE40B7jcHRD6/PEdiXOnhf06Q14DrJhkc0B2ys7Or
b4/Hvf3Xdj8bxLwEunppx+Yi/rQAsdCw3F1xIkJLoOdhuIy6xJAvIj2F08RJIDy+trgR+5iJr0fM
FH0G6piEfuHEgnmkvMH8hLqEg9feIQeQu0UkFzIMe/Lz8u3yePn469X6/AS8a++3yjzlF7osUEu/
qY4eMZdkLK5n8wlQEuIgOxeH/Hpe0l/gO/9/lqRogQgCSJACKeUWARiSKMZz+P3v6j9RtVQrUPs9
P2cT8tDwxQvgjk+snLbGTZnTexTJcseotgFLDXRFQOEdPdP8vbbh/NlBEPuPxCnV23eghGIyIxET
CJihBn9E/PBjGQP80VBYjBEF+JDiUPT4UeXcLE+UOiJxEvArlFogItmrNySiQxOqKkIMFYIIbOEW
cCIsoh+E9bCH4HbFwBvw84iWW2nswFBw+lU5ucMg5n2iepoT/5hooHj3gOdOBv9KKevUw99P1T5p
YHwBNoh1O6GgtaFEKd45Dq+KWRCnLEEDpZC7a3wmwRsfhv0CxowmtJWmVjzT4RNf1WFxgWggmuI0
cvNRxgGUBYRJEQIEBI/hyFsQ9aKgSi5/Du6/syxHYcaqjS5sOD++LIBBIQQkfbYUy15BndVDdAQS
g3qwDWAAFRCCJwgUgG+CXq4EBSGlGuykNvziOXLePK/MxtZ7gkRjrfaRpC+qxkxGsayxS/Q134rJ
oyppCRcq+3ljlY4YdtdYXNM7JhUbfKOSqKNisC7GSIsFbEC7cUxt3TS+uMvbOmkIwThKCyrFGCSI
Q11JgBaaLJvxIZMZGIhgb22lKJRWFq0hHFQwA13/UWF0PEDgBrzoyQhGJCSIkgpFiSQGQARFVYow
RWLBiEDNzCkbhw0osLpjyJ/CEGmqIqEgg3FyxUPs8AcQ6wryAWbQjImJCohUmvZ3gaDY2YEQElpI
RCkTXB4OzDoUHZsN8g4rygCyBIAdrMEdXWKavJU6RkY19wnqVT6w5m3KtT69xqKode4TqOJActeP
eqNlFYRIRhCIwIIAMCgqhWiBoDrnBJGE91iijyGQXy9aTFD7cxYmgAjiaQ5zI9w+DgF8j7iPkfqH
Br2pmUiYTiPYb8JfJVkCRQLAZ/fgbsVYCyAlWwKIZIDGoCSJRAK1+3r5fcb+GPICRDlTRr70/tO0
+Y6vNWv3QaR61VVVV/GwFoH+JQta3qOP9ndA+n3vwh+A47z6/6gLk+clGPzno0Tjfo36Jz+rOTfg
OVyfPxe2xprMwZFsx4TBucngYYdOu27N+s4N62RxaXkNDnRvg5Eo7I0cHRzRwc8DPIa37wO0yXvA
UnHjw/gsLcnoJDk43NdcpxxyU788JyaQbO3Bzc7mEU58YM8JyYkN8zo8Ewxw0GDbRgUAnURBZUT2
lSzJWGRjgGKFzJIVgkgURD6vdtsbHIzOZx58TBuHE7XyxkQ+BFkiSCHv9Es0IiIjGB5/jQoP6DJ+
OWeH6UPfev2zzw1lxchTQQs0pxHQYeIMHr50U6zrz/Qp3d4ejAPDzjkW8iayBIPn8FodvFd9UWXT
sDYqH7eiSDAyDI4d87bTYwCki0RiRMLbTCotNt7fOq79tMLBmApzcp6FF0sJtr1aiqB7stf0GcKG
ZptiNlDINg/aINBAdOPP4AJFDVeZ/DjR4FUxOFghygnKfggDHzCIRBX08lQKU1UMWBQB60lhTggH
lDqFwvMl9G+N0RmuAd8Xz+b8bWkYwjCJhVAmEQKYkMJXuU4ufyALmcF8ouhcQtrHQD37CIfMDEcB
4E6wkE9ZameRy4+WbC4ifLJtPZR+AibFMgfAGg8hz+glK3AiStApdkgUYExAhoPM4GBPmB+YBSbe
qsTntR8KA6I5IDgIhSqYor9bBQU1OHigdcuntFCg+ioaLvNlfWGsHHQ4xJELIJFQb3grcwA4cngB
6JU1vlqQgOg4EnAUFQSx252QbDESg1fmdvc2ChGoDvgJoMSAhAreEVQbyjslnAuDgBzVnxFMgA+o
OP3h11CUCPiaIlwTgUgeKoYgcBC0eaUrinS7FPiVD8FQ3eLqPLQMgQ8vCtNwYXoUWwoJaf60MVRR
JvaWgX2EDbBknzIyaH2/+Kr/EURBGiucSG5ZuSIGREDziQDYfOpQ6B6hTqH5VQ1Vw3SQgyOoiFRg
BpEFEntpFB8DsMBkUmNkEK9OSrH7XKwXoUql+menkNP3wP1mClxYhulE7cDeiDIqkihrAKghIkYr
IjIXXFbndspgwhAmwrMfqIMZECRQikFgjASRBFGIoyIoxFGQRCRYCMgAskkCKAyJCAqxCJ1hUAgR
VB3fuFxXGycHb9dg4bdUxg3PZhuAQsYQoNwgWRMTS98IZ1BhJEgxhkYc4/sEQxwg9SqVrEeMgvBv
Jr6lDoMDZlRrOMeoFc/R6fa8Y/oHk7PrMzcs95UVg9gQLEieXEx7vsw/JVBKoTIDM/WVf3Q6gwMb
9hdX1oeCuq0mBkwPnDVCET4/k0AoDhR+U6S94sTsYog0F69d/yIjYcNC1BI6kcLlxoMMkSyfxBxG
0kCQirCKBGEEOSJHCfMnnF7PAe0hv5mHtRT1i9QorwmO0OP4D3xrIXzeY5CEOwzYHb0JDjBOsQCP
Dqfd9KrBH+UWWtp6kPiwyKojwBIchhOuH7h0NeBfiAeyBPVdGMhP2yoQQkRiHeBlliilJ8njQ0BB
65uiiEgIsLhRSJA4+BJMfgpzFeQXRc9gQ0t04v3AQUphVHKFA2VVDtDPCcBJRR6RppmUahZkGJeP
hcLbrpeWEujyWYoFFJQWQYcy0A0SVskGSCwRPFKiFUUJmNJc/jyRy3fg+VVTzMbGYmU0+CoNwxhk
L8K+EEAH3z7GJEWEkEDcnoYpq8Usqbm74gIxQ5Mx9HNgaSahhuw7ZFVfittW2rZaq2yqqqtpXR9f
s51tzATt4sJOd/ft1Tl4sMzVVDUjCcSqk6TA3oY2zbZ2JbLyYaHo8d8cjNXe0mF7gqhYUGBAeHYP
YVzb5RzkmxDWKFiSSeUnXHdeofYswI4E3YXroQuXzIZxbY0qrDmApVqSMR0qghKoVEp0wRxJ5OoG
pMXW9uqsA/v0YyfUE2CFAxKDOZQnmGHAKMEgmGGwJvoMqAe9Nhf0C8Be+73kIQiRBGSetliCkFUl
GSpAiSMAgsFYiLEBILARqCm9NQzSBR5sAf62PRufZiLvIO9u7QDpeoGJIj1/vBfGSdggfCCqMFVk
RH3iURVFCHzoda93D4+V0OAUU4xAoNfMmMKjTtoWhjCSDCAkv/c1sMDTIP5QqO3Eqy+lMVDxpcje
ggNnPm6UHAQ8j2skh3RklwH64o9qGJAZOZp6lGdQ8wIfbLegf3UWygtJzM2IdK8sRFsV22HK929w
Ru+6hQhX5lm2Al/JSGuRmgjrnbi54c7jGA+UFYEFgAe1FPAptUPm+tV5R98QY7BEMRPsET8ADUQN
hxIeb8UbPTzv6+hyJ5pAkQkIYDIIwaAUdGk8oIXFNdFD4/lpzFi64fUFkLRGoR2urzY45G4vHeoU
cW/OPaoX4T5AD7UB2LuV9inBmqD4aaunsU4DvETDAnc/aKRAYgFetc1DuXzlv1VnTCYfuQxZEDSi
nZQTHAxGMQYiCHyso8EzoUA1qEnZ/OnqMzMTBXI5UcVbQ1rqjZLBIETVGgF3KoWFoF/sRgLhDhBZ
LxRoIpGFURYKbAUBdVYr5YqUjEjEhGAQdFEE3wOnfA+1IixDmVUfbzKGDgPDGBsgWdKjYgtuXrRS
grz7OXh3AHGu4J2TYsE6A6owjCPGmoxQvLRwhFR4AZRPM94DMB3AhqpGrc5QzemK5JMIZYbm8nGq
CbBFJQpJYG0LZgyBvcMoTcZJohA8e6CKkWRYqkm+geCpUlOfKBoFOLKJyhWc5DQIthLCksFbLjrx
TLjMQ8mAYnXUCQGkTUmuZTKBMs1gG6UfydiKnGqBF9AiEfOgPaG8DMOHCjHYM493XQDK/0ju3Qq1
/Z/LmtHV3s7hZaNlLZfVmZC0lAQkPAWIJQEAIhNWihXAWUNQhA5ssEg7uxcMwv+85ZfZqsUXIxqs
SZ9EFYwHSLKLj+QHGKUpcYQCQiR0KqQTrRxEsvR7bB/MYi4YAFUaoIxVXqJuQRLbXHP4vFT7Mtxh
21ejyWAx+WSWCBUfiw5KRUAHrP3CJ4bGd8SFnsKalMgRCoQCyYAF7MGA1auzh6bAJ41XaImSPpQH
QrQxSAllbAlgo1Nnw+Ou3ooL4UEYpFgQkBPzr+d9RTAgWEOQo6iKcA0Q7XSDCAdrsoUsPJA8H1DE
li2w2BQAduGEgwnMglkQpEpAYBQiWhSyRJKAliUiUgwKQSgykSwQoCWCFBKNgliWCUaCWtwuoBjR
QgUBy5gPJR1VSioj2KxKUC21sSsgx9sKWAIEkQNx6QbKZ2B0YMmDsHRGj7OkgVhm2LGxqzs2xS0N
v56fdo2gbBLAUaWv6fXZykQEE+5mIGPFCrRL/R/Zge0iHQqEOlmHoCcF7G5HpF2i+C+kkpH9iA3F
tcnQFrFpY/8WuiHT0i17U9Ym/vIXih0ob2lgJmnHwzsuad7pQszAeMBSChSG9dhCAB/CF0JdlQbF
hsoh4vsEx7fBi31UMd7VVBEN58wl6BydmXyuCBlKIIblAtFZaB1TITBIib2gG2qYI7FOtmoMBE0J
oijS2Jma42hix0lN8ijvoC6QTC+BYuxkX04WwuhCJTjKUIkbaYKZcUEF3OlJjNOwyOFP/0EmJSEG
GAIyKREgCRFCJJARcEs4KWQBVKMCvxSQpIkiSIiAMiMWSIvJJCokSRhGfjADH5B+BAc1fyDEcoqE
CEkkiQcWmcabDGbikBqAiSKJb1nrPwttiGJgDskGJgtxwAwFFaGCCQIoSChIThKArVcjLJEWz3lS
IYI5HEnHbFFQN3AWhsojPEwAC4WRIJUCgjNEEmzMEVQNaote+IaYMgTWwKoNImAsYiNxbV6xgNAA
jywUDBxGCoYiIZAjAC6MSAXlMGpLJCwbICwHoESYIHfbIHWFFuCMyAv741SPr69oZkC4ULqlZqmo
pssqLYkkIIsIiCN1A0t8bYBQuyRVP50DgIWa2i2sZXCuCjguG0GxhZD/WRSxWzyXsGHehiU+8gPD
+cB5hy9YPURkCQP1n6s/FZA/C6FiFPppS0taHrXsaKaaA9ELEGiEQgBxGxE5mqRkEBfsiA96Nugg
VZgtCRKAtBQk7gQJHAwCXuKSYqIMkUJNrWJRGoqGBFEKCKWiUrgNKJUUaYwkggVAEYCLEAjcgkY0
kwiELNoJIUQokAKiMAkyrIKkIMIMZCKwUijEgjgGCS6+CkAjAVYiQALpBEuNUBqkokBkgRh2nRgj
D0Fge9BEKmADgImQ0BmB1ghw48ROPbCQuWIkfck2HWyKbHmBmobpgAm20UobR3GTYOzt7/uP8h7/
EFKg3aMMyD97B/K7h5sRn2naQEaIbbWKCwwcSdj5tjnNjsnw9Vh51YK9oI5e4+X0HVbogZkAXcQE
eIkUFvj9XnpddYMiz2qH7VDQTpvF5w0AG0eCbXXueG7uKK503nEeGSHp/cjW8gZty3GQJJJJ/Psa
qA8Aj/ZQDY7leAFPIOQcASyfainl+/6wIEFiyEn20kKpEIqRggxGMRZCKoEChWMgDiIPnTEDx5Kg
UB0sYpgXf7xBCMkUIkGAMkCZzggiiArBAggyIkIgkYCMYCjGBESSRGQBBkFgoIisYhIgSIEEUEFg
kgDEYgqpGDGRBJBGAAiAgZvRFdxqoOqj417+Xkt0o4Av8wIv4gxFIEWEZBkGEQSSQkDU7VQEGRhE
QSByKUIWyeDkLTDApBIwWUYVgCpRawYHY5zBqk8Wg+2wgWKJJBRSKipILEh1F1MkLEgIyIqKoJId
ZSRjUEBFSQYICwEIQEJEYsIwFSAEQEkFDpVpRaAihCIsXe1Ak0HmEOXAcH+AuGABxCPlTcTKoQ5o
xNgacrA3yqpAPOwAqpCGmpbVFgiKSIQYQSDICMUBlKFxBXAWAJLlGrELEgTPxDqfuUAIoWRISKjm
B1GGg34jJyQIrtRzKlPq7CR6bwR1aAogB3FKJESDIoznQsU7QZswDTCn9Ps31FkScXZwmUhIBSBE
dyAcYNKpmQRrMCgBZJH9RoCrMEKJsTZ80Hsuoe/8zoORp8D/KMrBB2IB6AiQiQgKRiNAjtB53nDj
HrXiDbWgN1EN1nhTIgUh/UIFgOU5S1Cmjxo8IOAwvR8Rq333ecf6YB2qK+VOqxyUHgLK9KAfxAiO
1C/6HYkCCgeg+nwHjA4hAnI7ih1T1hn4ZQ/RnyGtVF+cP8tA+IXulAqm9D1lNow2roWj6jzVHTQE
I0caGfNPnuZH9ppSBbI0KS7djG7wdvlOA8J5PWVHqL836PkCB5zu8LMIX0UAp4EPWw+h/mK0kTOj
IclijgccjIARuj60Oa6dDykeUXMRDDBPqGqqEMSiiV4jHRvHZVBmQQuRhUKHhSIn5p2q4nYc6hkL
ZsLARYKDQVwRORh9v3YAKVrDfCE1GAAUopzoDbYWW64eGtKKqxgput1wsCZeWGSofSgOnvG8yJvG
gf/TF4v5rF9uvGgwBXADvAsjSd7A+4aULKpaC4IwWn6kERS1/UkT2kB7vfQRyiBchBAdhEPC0BIB
IBaPjCqIB2gU1lHRM027McEqVYqqqqqqqIqqyCvUwpwSWSk/pGVjIs5sCF+KeHZRgGrjZCQcgCUg
wkigJFdhSgFRFmuSWfrzWCRNi0RVYqqqqIqqAVEpkghOgcwOtwC/M/ECPliqewRDynzeQMBAP2D5
e64H61WYAK4EcNCnfUNUh0F7CD6Ko6YdYAUSQe3/7RURzVCFTzC81G9ykCf2PwA2ttNETIfwVcw2
M6DrAUfUekMlEOoCb2AEV5AsJzP4qp/cr2nyopcNB+NRMH+IidgPlQ63aJoKQ5QDrPeCwQgHekXx
/K18Ya2tEU31JyioIowFGCQD3lkoyDCMYIogsnvZVQWIDEZEikWIArBFjCONGfGyS6sowHTJUUGC
qxCKrEUVIPkZgojJDIeyMmMLBBgS1OjZTqYryI/oPhVTFAPWBg0DGMiRgxgZhAgBEdCKWFDzMB+6
wt9FC5C3KFyFLlNakRIkgRYMxXZVWVBgqHuY/W/fG6lKh1pKPoFwBSxMjwEUtFIrBCQBHAHIUXzH
MRO7XNASA5clazVjeqzAJQwc2LlTIsIluCTDDmSBhCFjCJBB2Tu+ohW02RnxX0wYcnZB3lVqKc6h
8Tyapw7z/mi9iA0bhPUgF0W4BcEuFkWyA2LCTgacN1w/pJgioGAgBIKDzH9EAvuQ3bY22pDz74fl
pdA7z+gUc4P2f03sJ54feL8haEANgOgc/SIhOqWgCHbAK+5OdVEvJAkjADsg34/EZhbP4C6hLodo
d5AgjtlvsqUEzrQv6bjfM/KZPPlityAGJSRVRpoQ+lGKAlhNUAhQMRGnUyEzc4bxIRgkUhFQ4PwI
LMEQHrCIgpsKBpEFuZkwoMJdAs9ebbM7amqS/g9hWt0JGwkfhPgoECB3UFrKu1i9hVuEAFOggXKs
0n8loLQooaYAKWCAUsQC6wEoCWiozyAocl6IEgAQvcH4CkMJuoackDASRKKONIKGBGQSSoLTUQFI
kVCcwIe0wtMfVzViGJwQx7nbJrLiUg3mmc+mZNgTnyOyGHCHURYmT7wpwEn9FYbGjRoF6Fx3G1cM
WsWfnlN6OItdy6EM3+dx4YUQlnZjEi94P0LuDbDyXqFqZvyC7m7G7I2rdbVJVB15V3IMtgxy8i40
hTlotbUZVKNZg+vpTcduXazjsdk4/c9rgXKptVj8KjCq7GlNWps6vgddGzky6HejkzIweq2RptAw
SGTJZkJMBSoYd0rAwsNzmZvyCTbguENJCaq2rILEWXxDUsmEuYY2egIIWAgUKGSeAQTaOAXHEzxz
VNaDwtxJlTvWbo+yeYtj4nK4kqaFWRMOgdAzZRAZ3wpwwMDgADnAw08AYEzBnr2Z2smkBELlGJjg
xqXU69rElrrm7tBmsRx0YyGSQhkkJkkG2NsbY2xTDBDEQRBsUUBBDj0oCyRFuoQDFyVC+CuKbluL
+G+3H8wuyrWnpTVHVvp3eL6CEhCKoqSIEAbEI28546fGfzPLuGMv5DzQJf6gwAkFFYEwWpbbUQFz
sR7qRr1KHc7ufUUepRVQgd4/gpcU0ggHZABDQBRg+XeyAKCkRgKAoCiyAsIqkXoqDrzFgTIu3jiC
hazyFnh2+YOjG+FU0VmR+YEIBH71VQjgqGwuhSCCcC55AS5HpQLUEAN9mhSC36hotEWRQTogOAge
LuNnxgEgSEgkgyAoCgpFFixR5kLCB8SHQ5hPbDQDeHTNPoXe0B60e77rWLWtYtalDzsjIzzvhCJs
xVP0/kv/amwcVQ94+2ecJNmKRQlKhBaPeGdsFFTbtDCBKTqbAqVKEdF4mgwkmk2KhsXQhhGDEipk
2BChhSyUZIIu4ghB2WaKBjAwSCESWSm6uDRSAsFGD4dOYO8PyJhP3BlGCd5jcWDGYKrmUFUQcLhm
Mo4Zhcctgqk9TFiKqoLFWIxYMGKCwBH+QZL5+rzn1duJAwOxUG1jzPNHdtri84wUYacdxYeChHyE
38abQ7gTuEQ/ah92OfHv6B7i557XD4aDbPHC9zOylkyRw8C5Zhl1LY4l8hUxCKiOwKryMaU6KaRY
MsA4CBaggqMUti4REtUZZb66p5UG1KmTLBS6rBnCCJsqjAMf31NJy2EdssgM9GI0gaUKi1bdsyqA
1Cpr8CYvfLNpsflUQcUvSb7YdCxwc1tlgfdyR5IYhGdmxIqHhNLVZz2JOs4KkJeDkwwI6YrTvaal
2woJ2n8dYM8bg7ogbCCiAymcVMdCFCwgiRQmENDNYEoyINwruspHPe4MqV5gyHC6dL0Rrk9w2Y8w
NbJ3l5EmXa4SWSjEZEElOG+UQshLbWJP2q/mfgD2ZbYu8gJrarmtcbVL2LNECspNZMO2aORI4EWA
s3KCGi+hIGpDQpMBoFoYRg2QQIOUtiQWkHiUIi4FmrrQBcisIqySDcEWCpdiREvSgS7cVgpHqgDk
xAG477u4d+6AHUILc1OZqCeJZ/GJPGdv46xzLjmWQPKLFjwTn0pDz4IsQRHqnIDZxDFUUgODIpIQ
tGok50GMsVGmWpt397ZtHMWbnJMyKMa6dblsCrDoQqxWNJ3Xq+MVPfpZUJfY9TJ0pOWVG0OkP2Jv
cScolBDzszNLWKqnZYhaBeePZRkRduAtBYzNHG4reCsRhehkECRURSMiwYCIwZFBBIUOwNd3InDy
mKCAa8IBCQIqaOBqSff3Vaw2GMD6E/XMm2betUVs0ZCyAeGG4mENAMQJBKITKWHcbgBShAHUXUE2
A2+xIqdsozIpkKDsgSEikgrIOyJz8zawXvS1H4Zl04B7g6o8cg12BguuQGkXhDlwEsglyEGMUIAa
XM/MC+K+BGdOzgp6m9IvVBhhuJZFGCIIyJ1WrBrS0oyKBy1RyQsggwYgk0liqUSwhEUisVihMCqR
MRUGgqwIY+0RDC81DnGo1BDsg1LRhAMzBK4AB3PrDMSPAYadr6DUA1IC8kkiwsQFIIJBBikjIyMA
Ug01QTIBGhhFDUWwMBegdDoGL0RLoKaAHCIbBOTFHoTHuAUOK2SJ/0ulKZl+k2HdfjaygFvOkdAL
wqNLrB7wfr1e+nIzf17bOt1IgqyxDQcBpIiYjduMgKYqSz0uKDRdihbEvkZminZFC0UBjmoHCu42
G88pFJEVYRigSXwfO/JGd5VPlmXDApSzBDJhQtopEtEsoVsGAiUs9JhdH8RIf5nyw+YGV6UB2Kve
j6Q9KRISJBSSCkRHxsThOX7gjIsGQInaAAZgeE0+dHoAtFRkD8e5BGh4vi1CFkRbayhNh3BtfpIc
AbtzkJ/jzI+5cnpPQTzdZYghp5jjOEEPj51EE8O5gu0OZfzggw8+SB0qHOaNMpDwov8IsjzguEFC
WpUG005ZofrQIOqoFhHxI/f7MsEM8H061L8yICcCJhEG8UQqFPypayrcjgUp0Za50GoCj6DCAfVD
F96odQ9oaGZoKDnLMLBA7iynEN0yYKv9Fi9wOvZbR1GJp7VMQeI+g9/pD/K4P1kjDyE7g0ElHBQU
RgxiZAn10AKH2QBQ+QT3I7QP7keYG6qe0cTgYcDDihISEwGgACT4gnVD5X6UA5ucgWkCT80iAaGW
R84epGghuShz3zN5hE/x3wuwA7ACCwGADFGKKZId5ilBoiLF3kglUVsLVsgeoSfQgQiDDTy8Ca/3
Kh9x2OUnojmGWQBxgAl1DQ06J3CfqQ8i7ID9o4ZhR+cI0knmD0AhT1+N3+dDBIJH9PBLhGAwWCyI
DhCiCfufsQG2O+j9KKe0RMR5dhk6ngoLiSHOahLFETYRHzmxAOcEPkV5lDB1+Jo9kgO2uiG6fQ+P
i9l+dIRAj4lyoMDkaO2K1HrjkMc0MeSwp7ynkE8yvQgLzVRMLmVZ8hH7usyLkgsuZdf9LLNofs+X
siRZGG9DAgB4sGEREQSC+IYJBlGc/YpYs9VPSB5xS5vAGMEAjGAgD2uhGEK0IAEgSKvwqETUGAbh
Z2dgiFUEH7lP+nE7+5D0qGKqj+8CKHR5R6DyL/BQ8h37i5Z6L8YODc8JFkgWtnydeOK4Z0Zi1C9Z
y+y92+ITmq8VMr0r0l7HspEET6bxFzcMGMyqI5gxaoz4WWXOg5RM69rCUsJRPlClB5aKFORhRtiY
E9Pq2NTrx1DJSggtqw2EbFr70FyaeF3RIodybV7FOiH4PAdwO/SE4oi7WECIkA8SIoFAWwKFL6q4
uKofRQPFseZDu7CQ5KFUamAggO4bqa+S/ChrWwiOthE5VekPbQfghp0CIdXQhLE3c9l4+kEpHRI5
AeZQ+VXiFLmCdBnOvoeR8QlKI9U51QoXah2g2VD50B8GnyKVvb/vQogQnygJ9gqsVEJpkBEZC5P5
jUyQEYJpEu239ybQJ+CIwjGQEiRQP8PwpnCCCzP/V0/j682L4/lNxaP1qfj12/3b5bTR1BdBGEWC
QsdhGYVBKqsapcB1sl0Ngo770ZUkBRUbSg6ppkH07myUGyerY/ClN+/9/6dRP0/Vm9FIwRn+T+5e
y6pdIzbcfzosvGMPd1aJwwOHl//Ah8wQ//i7kinChIJcw/so
Follow ups
