diff --git a/README.md b/README.md
index 27cac92..88cfbef 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,5 @@
-#Introduction
-Welcome to FPAKC (Finite element PArticle Kinetic Code), a modern object oriented Fortran open-source code for particle simulations of plasma and gases. This code works by simulating charged and neutral particles, following their trajectories, collisions and boundary conditions imposed by the user.
+# Introduction
+Welcome to **fpakc** (Finite element PArticle Kinetic Code), a modern object oriented Fortran open-source code for particle simulations of plasma and gases. This code works by simulating charged and neutral particles, following their trajectories, collisions and boundary conditions imposed by the user.
One of our aims is to make a code easy to maintain as well as easy to use by a variety of reserchers and students.
@@ -7,14 +7,14 @@ This code is currenlty in very early steps of development.
The code aims to be easy to maintain and easy to use, allowing its application from complex problems to easy examples that can be used, for example, as teaching exercises.
-Parallelization techniques such as OpenMP, MPI will be used to distribute the cpu load. We aim to make FPAKC GPU compatible in the future.
+Parallelization techniques such as OpenMP, MPI will be used to distribute the cpu load. We aim to make fpakc GPU compatible in the future.
-FPAKC makes use of finite elements to generate meshes in complex geometries. Particle properties are deposited in the nodes and cells of the mesh. The electromagnetic field, with the boundary conditions imposed by the user, is solved also in this mesh.
+The codefpakc makes use of finite elements to generate meshes in complex geometries. Particle properties are deposited in the nodes and cells of the mesh. The electromagnetic field, with the boundary conditions imposed by the user, is solved also in this mesh.
-#User Manual
+# User Manual
You will find the user manual in the *doc* folder.
-#Installation
+# Installation
To install the software ...
diff --git a/doc/logos/PPartiC.eps b/doc/logos/PPartiC.eps
deleted file mode 100644
index 5517180..0000000
--- a/doc/logos/PPartiC.eps
+++ /dev/null
@@ -1,177 +0,0 @@
-%!PS-Adobe-3.0 EPSF-3.0
-%%Creator: cairo 1.16.0 (https://cairographics.org)
-%%CreationDate: Wed Oct 07 12:01:13 2020
-%%Pages: 1
-%%DocumentData: Clean7Bit
-%%LanguageLevel: 2
-%%BoundingBox: 6 7 205 57
-%%EndComments
-%%BeginProlog
-50 dict begin
-/q { gsave } bind def
-/Q { grestore } bind def
-/cm { 6 array astore concat } bind def
-/w { setlinewidth } bind def
-/J { setlinecap } bind def
-/j { setlinejoin } bind def
-/M { setmiterlimit } bind def
-/d { setdash } bind def
-/m { moveto } bind def
-/l { lineto } bind def
-/c { curveto } bind def
-/h { closepath } bind def
-/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
- 0 exch rlineto 0 rlineto closepath } bind def
-/S { stroke } bind def
-/f { fill } bind def
-/f* { eofill } bind def
-/n { newpath } bind def
-/W { clip } bind def
-/W* { eoclip } bind def
-/BT { } bind def
-/ET { } bind def
-/BDC { mark 3 1 roll /BDC pdfmark } bind def
-/EMC { mark /EMC pdfmark } bind def
-/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
-/Tj { show currentpoint cairo_store_point } bind def
-/TJ {
- {
- dup
- type /stringtype eq
- { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
- } forall
- currentpoint cairo_store_point
-} bind def
-/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
- cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
-/Tf { pop /cairo_font exch def /cairo_font_matrix where
- { pop cairo_selectfont } if } bind def
-/Td { matrix translate cairo_font_matrix matrix concatmatrix dup
- /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
- /cairo_font where { pop cairo_selectfont } if } bind def
-/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
- cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
-/g { setgray } bind def
-/rg { setrgbcolor } bind def
-/d1 { setcachedevice } bind def
-/cairo_data_source {
- CairoDataIndex CairoData length lt
- { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def }
- { () } ifelse
-} def
-/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def
-/cairo_image { image cairo_flush_ascii85_file } def
-/cairo_imagemask { imagemask cairo_flush_ascii85_file } def
-%%EndProlog
-%%BeginSetup
-%%EndSetup
-%%Page: 1 1
-%%BeginPageSetup
-%%PageBoundingBox: 6 7 205 57
-%%EndPageSetup
-q 6 7 199 50 rectclip
-1 0 0 -1 0 71 cm q
-0 g
-49.766 27.305 m 49.766 29.898 49.078 31.813 47.703 33.039 c 46.328 34.27
- 44.297 34.883 41.609 34.883 c 36.625 34.883 l 36.625 47.945 l 32.438 47.945
- l 32.438 14.352 l 41.516 14.352 l 44.234 14.352 46.285 14.953 47.672 16.148
- c 49.066 17.348 49.766 19.242 49.766 21.836 c h
-45.594 21.398 m 45.594 20.18 45.313 19.301 44.75 18.758 c 44.195 18.219
- 43.313 17.945 42.094 17.945 c 36.625 17.945 l 36.625 31.289 l 42.094 31.289
- l 43.313 31.289 44.195 31.008 44.75 30.445 c 45.313 29.883 45.594 29 45.594
- 27.789 c h
-45.594 21.398 m f
-75.641 27.305 m 75.641 29.898 74.953 31.813 73.578 33.039 c 72.203 34.27
- 70.172 34.883 67.484 34.883 c 62.5 34.883 l 62.5 47.945 l 58.313 47.945
- l 58.313 14.352 l 67.391 14.352 l 70.109 14.352 72.16 14.953 73.547 16.148
- c 74.941 17.348 75.641 19.242 75.641 21.836 c h
-71.469 21.398 m 71.469 20.18 71.188 19.301 70.625 18.758 c 70.07 18.219
- 69.188 17.945 67.969 17.945 c 62.5 17.945 l 62.5 31.289 l 67.969 31.289
- l 69.188 31.289 70.07 31.008 70.625 30.445 c 71.188 29.883 71.469 29 71.469
- 27.789 c h
-71.469 21.398 m f
-99.891 47.945 m 99.473 47.945 99.016 47.82 98.516 47.57 c 98.023 47.313
- 97.602 46.957 97.25 46.508 c 96.508 47.469 95.535 47.945 94.328 47.945
-c 90.047 47.945 l 87.617 47.945 85.852 47.461 84.75 46.492 c 83.645 45.516
- 83.094 43.836 83.094 41.461 c 83.094 40.508 l 83.094 35.551 85.426 33.07
- 90.094 33.07 c 96.047 33.07 l 96.047 30.664 l 96.047 29.676 95.742 28.91
- 95.141 28.367 c 94.535 27.816 93.641 27.539 92.453 27.539 c 85.016 27.539
- l 85.016 23.945 l 91.969 23.945 l 94.75 23.945 96.801 24.555 98.125 25.773
- c 99.457 26.984 100.125 28.938 100.125 31.633 c 100.125 43.148 l 100.125
- 43.629 100.305 43.988 100.672 44.227 c 101.047 44.469 101.664 44.586 102.531
- 44.586 c 102.531 47.945 l h
-93.844 44.352 m 94.801 44.352 95.406 44.129 95.656 43.68 c 95.914 43.234
- 96.047 42.656 96.047 41.945 c 96.047 36.664 l 90.047 36.664 l 89.211 36.664
- 88.523 36.941 87.984 37.492 c 87.441 38.035 87.172 38.719 87.172 39.539
- c 87.172 41.945 l 87.172 42.781 87.363 43.391 87.75 43.773 c 88.133 44.16
- 88.742 44.352 89.578 44.352 c h
-93.844 44.352 m f
-114.047 44.352 m 114.047 27.539 l 110.547 27.539 l 110.547 23.945 l 117.891
- 23.945 l 117.891 25.242 l 118.953 24.379 120.234 23.945 121.734 23.945
-c 126.672 23.945 l 126.672 27.789 l 121.25 27.789 l 120.227 27.789 119.461
- 28.117 118.953 28.773 c 118.441 29.43 118.172 30.094 118.141 30.758 c 118.141
- 44.352 l 123.359 44.352 l 123.359 47.945 l 110.547 47.945 l 110.547 44.352
- l h
-114.047 44.352 m f
-146.359 47.945 m 144.723 47.945 143.473 47.422 142.609 46.367 c 141.754
- 45.305 141.328 43.91 141.328 42.18 c 141.328 27.539 l 136.766 27.539 l
-136.766 23.945 l 141.328 23.945 l 141.328 18.18 l 145.406 18.18 l 145.406
- 23.945 l 152.406 23.945 l 152.406 27.539 l 145.406 27.539 l 145.406 42.43
- l 145.406 43.129 145.578 43.625 145.922 43.914 c 146.273 44.207 146.82
-44.352 147.563 44.352 c 152.406 44.352 l 152.406 47.945 l h
-146.359 47.945 m f
-172.719 14.352 m 172.719 18.664 l 168.156 18.664 l 168.156 14.352 l h
-161.344 27.539 m 161.344 23.945 l 172.469 23.945 l 172.469 44.352 l 179.578
- 44.352 l 179.578 47.945 l 168.391 47.945 l 168.391 27.539 l h
-161.344 27.539 m f
-197.25 47.945 m 194.469 47.945 192.332 47.27 190.844 45.914 c 189.352 44.551
- 188.609 42.508 188.609 39.789 c 188.609 22.508 l 188.609 19.789 189.352
- 17.754 190.844 16.398 c 192.332 15.035 194.469 14.352 197.25 14.352 c 204.063
- 14.352 l 204.063 17.992 l 196.766 17.992 l 195.578 17.992 194.613 18.348
- 193.875 19.055 c 193.145 19.754 192.781 20.598 192.781 21.586 c 192.781
- 40.695 l 192.781 41.688 193.145 42.535 193.875 43.242 c 194.613 43.953
-195.578 44.305 196.766 44.305 c 204.063 44.305 l 204.063 47.945 l h
-197.25 47.945 m f
-0.831373 0 0 rg
-17.277 53.707 m 17.277 56.629 14.91 58.996 11.988 58.996 c 9.066 58.996
- 6.695 56.629 6.695 53.707 c 6.695 50.785 9.066 48.414 11.988 48.414 c 14.91
- 48.414 17.277 50.785 17.277 53.707 c f
-0 0 0.501961 rg
-23.906 61.84 m 23.906 62.969 22.992 63.883 21.863 63.883 c 20.734 63.883
- 19.82 62.969 19.82 61.84 c 19.82 60.715 20.734 59.801 21.863 59.801 c 22.992
- 59.801 23.906 60.715 23.906 61.84 c f
-0.596078 g
-0.751181 w
-0 J
-0 j
-[] 0.0 d
-4 M q 1 0 0 1 0 0 cm
-16.941 50.758 m 29.465 47.813 l S Q
-27.848 47.258 m 30.453 47.578 l 28.262 49.023 l 28.523 48.41 28.355 47.699
- 27.848 47.258 c h
-27.848 47.258 m f*
-0.137102 w
-1 j
-q -1 0.235294 -0.235294 -1 0 0 cm
--15.851 -50.987 m -18.248 -51.872 l -15.849 -52.753 l -16.234 -52.23 -16.233
- -51.519 -15.851 -50.987 c h
--15.851 -50.987 m S Q
-0.6 g
-0.751178 w
-0 j
-q 1 0 0 1 0 0 cm
-23.234 59.613 m 30.398 49.301 l S Q
-28.824 49.969 m 30.973 48.461 l 30.313 51 l 30.098 50.375 29.496 49.957
- 28.824 49.969 c h
-28.824 49.969 m f*
-0.115667 w
-1 j
-q -0.694805 1 -1 -0.694805 0 0 cm
-20.193 -42.855 m 18.17 -43.597 l 20.191 -44.342 l 19.87 -43.904 19.87 -43.302
- 20.193 -42.855 c h
-20.193 -42.855 m S Q
-Q Q
-showpage
-%%Trailer
-end
-%%EOF
diff --git a/doc/logos/PPartiC.svg b/doc/logos/PPartiC.svg
deleted file mode 100644
index 9966134..0000000
--- a/doc/logos/PPartiC.svg
+++ /dev/null
@@ -1,159 +0,0 @@
-
-
diff --git a/doc/logos/fpakc.eps b/doc/logos/fpakc.eps
new file mode 100644
index 0000000..667a7f0
--- /dev/null
+++ b/doc/logos/fpakc.eps
@@ -0,0 +1,294 @@
+%!PS-Adobe-3.0 EPSF-3.0
+%%Creator: cairo 1.16.0 (https://cairographics.org)
+%%CreationDate: Sat Nov 21 21:23:19 2020
+%%Pages: 1
+%%DocumentData: Clean7Bit
+%%LanguageLevel: 2
+%%BoundingBox: 8 15 191 78
+%%EndComments
+%%BeginProlog
+50 dict begin
+/q { gsave } bind def
+/Q { grestore } bind def
+/cm { 6 array astore concat } bind def
+/w { setlinewidth } bind def
+/J { setlinecap } bind def
+/j { setlinejoin } bind def
+/M { setmiterlimit } bind def
+/d { setdash } bind def
+/m { moveto } bind def
+/l { lineto } bind def
+/c { curveto } bind def
+/h { closepath } bind def
+/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
+ 0 exch rlineto 0 rlineto closepath } bind def
+/S { stroke } bind def
+/f { fill } bind def
+/f* { eofill } bind def
+/n { newpath } bind def
+/W { clip } bind def
+/W* { eoclip } bind def
+/BT { } bind def
+/ET { } bind def
+/BDC { mark 3 1 roll /BDC pdfmark } bind def
+/EMC { mark /EMC pdfmark } bind def
+/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
+/Tj { show currentpoint cairo_store_point } bind def
+/TJ {
+ {
+ dup
+ type /stringtype eq
+ { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
+ } forall
+ currentpoint cairo_store_point
+} bind def
+/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
+ cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
+/Tf { pop /cairo_font exch def /cairo_font_matrix where
+ { pop cairo_selectfont } if } bind def
+/Td { matrix translate cairo_font_matrix matrix concatmatrix dup
+ /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
+ /cairo_font where { pop cairo_selectfont } if } bind def
+/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
+ cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
+/g { setgray } bind def
+/rg { setrgbcolor } bind def
+/d1 { setcachedevice } bind def
+/cairo_data_source {
+ CairoDataIndex CairoData length lt
+ { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def }
+ { () } ifelse
+} def
+/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def
+/cairo_image { image cairo_flush_ascii85_file } def
+/cairo_imagemask { imagemask cairo_flush_ascii85_file } def
+%%EndProlog
+%%BeginSetup
+%%BeginResource: font f-0-0
+%!FontType1-1.1 f-0-0 1.0
+11 dict begin
+/FontName /f-0-0 def
+/PaintType 0 def
+/FontType 1 def
+/FontMatrix [0.001 0 0 0.001 0 0] readonly def
+/FontBBox {30 -240 640 735 } readonly def
+/Encoding 256 array
+0 1 255 {1 index exch /.notdef put} for
+dup 97 /a put
+dup 99 /c put
+dup 102 /f put
+dup 107 /k put
+dup 112 /p put
+readonly def
+currentdict end
+currentfile eexec
+f983ef0097ece636fb4a96c74d26ab84185f6dfa4a16a7a1c27bbe3f1156aea698df336d20b467
+b10e7f33846656653c5ac6962759d3056cbdb3190bac614b984bf5a132dc418192443014ba63de
+800d392b6fea026574bb2535fd7bb5338f35bf15a88ea328fdaa49670c7852e3d060f3c5d6b07f
+2ef6d0f22646c5d18e19a2ae3ee120390f6dd96f76dcf1e127de5e9299077a00c17c0d71e36e5b
+9d5ec58fceda57739a6a4214d4b79d6c48d2784b60c320323c7acddddf34db833cac0cf109f799
+69d114a330d372e5c978a66acc84e3fe5557f6240856a013ffaa0199444e5c5036f775eba4a5c5
+8cde66cf604b9aca2178431127b8a1ff7ed633a65c04600af5f573483112251caad907dcd8c61e
+23500065b1568be79a17d2379b63d656e8f7715cdfdf357f0e30d9ab91d113c88231d875d60897
+1dee27eb5da34a09d2f1b1a2e7ab8be1036393cd6ae53eff0d77f452e9bf45eccc9e4836ae77c5
+b74262fa43bfccfd9a147a18dcdae6e50cc518129a64f669a5fae69c8082dec571d7d01a4e9f05
+6d3e52de0ea004acdbd1b3bf9b19aa17f839a75365a77b7275442a967093ffdf1694a72f9978a2
+304ac331d401b48e3e62a3a92cd39516dd480f088980d1ad8993a1f6fefb07e0d86b6f0293bb41
+68ac465726267cacb7516a0e910fe67c2dbfef06d8b64a9811506650d32fa182a0adcab8e2e21e
+ca6d0dc81959c25ea2d3f7ccec13e0cb4a7ef88e97c36e74fa13010220d6835ebdcbabdb507d84
+239e5483e8a8b7a52d6e1ea4ea1f5e6bef4534710c4055265aaa86fb445f3b2fc62cfdd9e283d8
+8bd083d09f0971cde00f2031b58b304d5f647f02aabf7ba9062c33979cd391f692c72ee179b7a8
+16f9c9e668d20021bd2b6a0f0114898729c6228be2895a696aaef0ebbcc842e64d5e72cc1d9b75
+44314028987a238f8fc4c18a0db3546c9ea42194b6bbdc45587e36d605fe2b7608d9292ddc0c9b
+be3e420b36fd52f9aef97f13533e101f34d4f882848f4845a7a824da815a710abe11a1ec8363c1
+06daa18dffb5a451af7bf3b20d79a63c94e305050ea893b7a61cc8cfdb2e8a593a073c2021b298
+40863c70742ab1734e1d6811cf1927832da10f562f6895575b50044179588ddd9ec4c413f68c3e
+3063f9594dc94115af4d9d4d6259c2ebb5afa796131772de3d297a8cc04f7f10398acc9142b1aa
+2da9741ad314918ff1553dafe4751b4c0efdc9d6a549acbf1b3d209f6ebe8f6561d627f37bbce1
+7213b92bf332c27718ca9f868f1724cda0774ea4c3a5a2ba99509eb9128c456e5526f234dc3adc
+37ac61ead9dafba1b5d58a9443ceb92474535cd3515e9ce357420b230fe927e81f06b2363c70aa
+b6e00858a44972ad3f8759069235bba0b8ae2c65a59fe3ee5642f88a8550a765907eb4f9432ac4
+9e896114d0bc969bc2c14acd9a50c31e2095133b6b4fc11a1136dbba4b515eaabf0cba23ffe795
+3532a1fca89780a841f3a5fe2514d31fd6d41fbcb5b8caadd53c1fa7b06506963f37006269d0ae
+fc1d5d6bd7f6788544e01a77bdf35aaacdd41dd4fc16237759516c60ee9b57e7b56606e0fb5a25
+9e6cb3b2e22d3ac4db73c228fda1b327cace7cbe25594ea2a9445efeda7826604e3daf18ca977a
+9a2788cddcea95b5460b648c12bebc3c39302a07481fe18ceed4c9e12ac7d51f8104bd589cf9f6
+68371d3aa8510b5ff08c972b84fd4e48c356f9098b9d130085f59bc7e748ac59e9715060d29b38
+d828e479d8a85e0aaaa7b4ab34b547643b9c3417305b29be240ba9514cd4656079f26f378a6af9
+f7c694b1332a14baa203faa5c6745eea7c9cc0b3fc13f722ff941857d141db5665b32ada8b7831
+06f5c01361bd4b161450546ad8812d364b394658c4d3e1e2cf1daba58315787fbe299a27a387f2
+1545b5d4bc3022e6a6803e971698504ebc44c9f7c02256a5ae01c9fb3194baf05fb7c9cd862427
+95ca9cfd57ba76d6697fc842ee19ff35ce8b9d3299a185a4d28537058b5ea471a818f521183ffa
+4b171f078577707cfbb3b2c6de4659b6dbacb952a0213b0586bd5879096ec53540ad07463e5b25
+1a38d715de6a0bb6f92d3e20e1b25be654fe7aad2e057428b68e2d6830f2ee0ca384fb862813ba
+d48832ec89991fa5e6bf9d3135378f0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+cleartomark
+%%EndResource
+%%EndSetup
+%%Page: 1 1
+%%BeginPageSetup
+%%PageBoundingBox: 8 15 191 78
+%%EndPageSetup
+q 8 15 183 63 rectclip
+1 0 0 -1 0 86 cm q
+0 g
+BT
+56.000126 0 0 -56.000126 52.419201 52.791732 Tm
+/f-0-0 1 Tf
+[(fpak)10(c)]TJ
+ET
+0.749999 w
+0 J
+0 j
+[] 0.0 d
+4 M q 1 0 0 1 0 0 cm
+8.996 27.238 m 20.621 9.359 l 34.656 9.363 l 34.621 14.746 l 25.941 18.195
+ l 17.332 33.629 l 17.41 55.852 l 9.316 55.785 l h
+8.996 27.238 m S Q
+q 1 0 0 1 0 0 cm
+21.156 69.262 m 21.156 38.387 l 40.305 38.387 l 40.305 53.027 l 27.457
+53.027 l 27.457 69.52 l h
+21.156 69.262 m S Q
+q 1 0 0 1 0 0 cm
+27.898 43.988 m 27.898 47.73 l 33.426 47.73 l 33.426 43.898 l h
+27.898 43.988 m S Q
+0.283465 w
+q 1 0 0 1 0 0 cm
+9.094 50.535 m 17.414 50.297 l S Q
+q 1 0 0 1 0 0 cm
+9.34 43.07 m 17.094 44.93 l S Q
+q 1 0 0 1 0 0 cm
+9.43 34.727 m 17.348 41.363 l S Q
+q 1 0 0 1 0 0 cm
+17.332 33.629 m 8.996 27.238 l S Q
+q 1 0 0 1 0 0 cm
+13.297 31.055 m 11.32 36.062 l 12.719 43.309 l 14.023 50.113 l 12.02 55.516
+ l S Q
+q 1 0 0 1 0 0 cm
+17.332 33.629 m 12.203 21.801 l S Q
+q 1 0 0 1 0 0 cm
+21.246 26.223 m 12.203 21.801 l S Q
+q 1 0 0 1 0 0 cm
+15.379 16.828 m 21.246 26.223 l S Q
+q 1 0 0 1 0 0 cm
+23.562 22.152 m 15.379 16.828 l S Q
+q 1 0 0 1 0 0 cm
+25.941 18.195 m 15.379 16.828 l S Q
+q 1 0 0 1 0 0 cm
+13.297 31.055 m 14.234 26.23 l S Q
+q 1 0 0 1 0 0 cm
+14.234 26.23 m 21.246 26.223 l S Q
+q 1 0 0 1 0 0 cm
+25.941 18.195 m 17.918 13.906 l S Q
+q 1 0 0 1 0 0 cm
+20.621 9.359 m 25.941 18.195 l S Q
+q 1 0 0 1 0 0 cm
+34.656 9.363 m 25.941 18.195 l S Q
+q 1 0 0 1 0 0 cm
+25.941 18.195 m 26.527 9.582 l S Q
+1 0 0 rg
+23.523 67.027 m 23.523 67.027 23.523 67.031 23.523 67.031 c 23.52 67.031
+ 23.52 67.027 23.52 67.027 c 23.52 67.027 23.52 67.027 23.523 67.027 c h
+23.523 67.027 m f
+0 g
+0.749999 w
+q 1 0 0 1 0 0 cm
+23.523 67.027 m 23.523 67.027 23.523 67.031 23.523 67.031 c 23.52 67.031
+ 23.52 67.027 23.52 67.027 c 23.52 67.027 23.52 67.027 23.523 67.027 c h
+23.523 67.027 m S Q
+1 0 0 rg
+24.66 63.297 m 24.66 63.297 24.656 63.297 24.656 63.297 c 24.656 63.293
+ 24.66 63.293 24.66 63.293 c 24.66 63.293 24.66 63.293 24.66 63.297 c h
+24.66 63.297 m f
+0 g
+q 1 0 0 1 0 0 cm
+24.66 63.297 m 24.66 63.297 24.656 63.297 24.656 63.297 c 24.656 63.293
+ 24.66 63.293 24.66 63.293 c 24.66 63.293 24.66 63.293 24.66 63.297 c h
+24.66 63.297 m S Q
+0 0 1 rg
+23.797 64.039 m 23.797 64.273 23.609 64.465 23.371 64.465 c 23.137 64.465
+ 22.949 64.273 22.949 64.039 c 22.949 63.805 23.137 63.613 23.371 63.613
+ c 23.609 63.613 23.797 63.805 23.797 64.039 c h
+23.797 64.039 m f
+1 0 0 rg
+37.844 42.566 m 37.844 42.945 37.535 43.254 37.156 43.254 c 36.773 43.254
+ 36.465 42.945 36.465 42.566 c 36.465 42.184 36.773 41.875 37.156 41.875
+ c 37.535 41.875 37.844 42.184 37.844 42.566 c h
+37.844 42.566 m f
+0 g
+0.141732 w
+q 1 0 0 1 0 0 cm
+35.977 42.188 m 31.02 39.777 l S Q
+31.531 40.027 m 31.91 39.895 l 30.895 39.715 l 31.66 40.406 l h
+31.531 40.027 m f*
+0.0679747 w
+q 1 0.486446 -0.486446 1 0 0 cm
+41.243 19.965 m 41.497 19.708 l 40.605 19.963 l 41.496 20.22 l h
+41.243 19.965 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+29.828 39.73 m 23.582 45.012 l S Q
+24.016 44.645 m 24.047 44.246 l 23.473 45.102 l 24.414 44.68 l h
+24.016 44.645 m f*
+0.0577315 w
+q 1 -0.845211 0.845211 1 0 0 cm
+-8.002 37.881 m -7.787 37.664 l -8.544 37.88 l -7.787 38.098 l h
+-8.002 37.881 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+23.82 46.355 m 25.027 51.145 l S Q
+24.887 50.594 m 24.543 50.387 l 25.059 51.281 l 25.094 50.25 l h
+24.887 50.594 m f*
+0.0732958 w
+q -0.252177 -1 1 -0.252177 0 0 cm
+-53.469 11.403 m -53.193 11.129 l -54.156 11.402 l -53.195 11.679 l h
+-53.469 11.403 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+23.453 63.156 m 26.16 58.465 l S Q
+25.875 58.953 m 25.98 59.34 l 26.23 58.34 l 25.488 59.059 l h
+25.875 58.953 m f*
+0.0654792 w
+q -0.576787 1 -1 -0.576787 0 0 cm
+33.038 -44.931 m 33.282 -45.177 l 32.424 -44.932 l 33.284 -44.686 l h
+33.038 -44.931 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+26.215 57.066 m 26.281 52.215 l S Q
+26.273 52.785 m 26.555 53.07 l 26.285 52.074 l 25.984 53.062 l h
+26.273 52.785 m f*
+0.0755832 w
+q -0.0137977 1 -1 -0.0137977 0 0 cm
+52.413 -26.997 m 52.694 -27.282 l 51.702 -26.999 l 52.694 -26.711 l h
+52.413 -26.997 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+26.5 50.707 m 34.664 48.957 l S Q
+34.109 49.074 m 33.895 49.41 l 34.805 48.926 l 33.773 48.859 l h
+34.109 49.074 m f*
+0.0739077 w
+q -1 0.214603 -0.214603 -1 0 0 cm
+-22.54 -53.911 m -22.266 -54.188 l -23.235 -53.912 l -22.263 -53.637 l
+h
+-22.54 -53.911 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+36.809 49.66 m 38.371 43.828 l S Q
+38.223 44.379 m 38.426 44.727 l 38.406 43.691 l 37.879 44.578 l h
+38.223 44.379 m f*
+0.0730145 w
+q -0.26796 1 -1 -0.26796 0 0 cm
+31.85 -46.757 m 32.123 -47.034 l 31.163 -46.757 l 32.122 -46.486 l h
+31.85 -46.757 m S Q
+Q Q
+showpage
+%%Trailer
+end
+%%EOF
diff --git a/doc/logos/fpakc.svg b/doc/logos/fpakc.svg
new file mode 100644
index 0000000..3412998
--- /dev/null
+++ b/doc/logos/fpakc.svg
@@ -0,0 +1,331 @@
+
+
diff --git a/makefile b/makefile
index 5b76fff..b7d7daa 100644
--- a/makefile
+++ b/makefile
@@ -18,7 +18,7 @@ VAR := ""
FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -g
#Output file
-OUTPUT = FPAKC
+OUTPUT = fpakc
export
diff --git a/runs/cylFlow/input.json b/runs/cylFlow/input.json
index 440fb58..e0d30dd 100644
--- a/runs/cylFlow/input.json
+++ b/runs/cylFlow/input.json
@@ -8,7 +8,7 @@
"geometry": {
"type": "2DCyl",
"meshType": "gmsh",
- "meshFile": "mesh.msh"
+ "meshFile": "meshTria.msh"
},
"species": [
{"name": "Argon", "type": "neutral", "mass": 6.633e-26, "weight": 5.0e7}
@@ -42,7 +42,7 @@
},
"parallel": {
"OpenMP":{
- "nThreads": 24
+ "nThreads": 8
}
}
}
diff --git a/src/FPAKC.f95 b/src/fpakc.f95
similarity index 97%
rename from src/FPAKC.f95
rename to src/fpakc.f95
index 0258a75..1482dbf 100644
--- a/src/FPAKC.f95
+++ b/src/fpakc.f95
@@ -1,5 +1,5 @@
! FPAKC main program
-PROGRAM FPAKC
+PROGRAM fpakc
USE moduleInput
USE moduleErrors
USE OMP_LIB
@@ -22,7 +22,7 @@ PROGRAM FPAKC
CALL getarg(1, arg1)
inputFile = TRIM(arg1)
!If no input file is declared, a critical error is called
- IF (inputFile == "") CALL criticalError("No input file provided", "PPartiC")
+ IF (inputFile == "") CALL criticalError("No input file provided", "fpakc")
!Reads the json configuration file
CALL readConfig(inputFile)
@@ -90,5 +90,5 @@ PROGRAM FPAKC
END DO
!$OMP END PARALLEL
-END PROGRAM FPAKC
+END PROGRAM fpakc
diff --git a/src/modules/moduleInput.f95 b/src/modules/moduleInput.f95
index 4f20b1f..7c24466 100644
--- a/src/modules/moduleInput.f95
+++ b/src/modules/moduleInput.f95
@@ -138,6 +138,9 @@ MODULE moduleInput
!Allocate the solver for charged pusher 2D Cylindrical
ALLOCATE(solver, source=solverCylCharged())
+ CASE DEFAULT
+ CALL criticalError('Solver ' // solverType // ' not found','readCase')
+
END SELECT
!Convert simulation time to number of iterations
diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95
index 69895db..7f3a4f3 100644
--- a/src/modules/moduleMeshCyl.f95
+++ b/src/modules/moduleMeshCyl.f95
@@ -10,9 +10,9 @@ MODULE moduleMeshCyl
REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /)
REAL(8), PARAMETER:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /)
- REAL(8), PARAMETER:: xi1Tria(1:3) = (/ 1.D0/6.D0, 2.D0/3.D0, 1.D0/6.D0 /)
- REAL(8), PARAMETER:: xi2Tria(1:3) = (/ 1.D0/6.D0, 1.D0/6.D0, 2.D0/3.D0 /)
- REAL(8), PARAMETER:: wTria(1:3) = (/ 1.D0/6.D0, 1.D0/6.D0, 1.D0/6.D0 /)
+ REAL(8), PARAMETER:: xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /)
+ REAL(8), PARAMETER:: xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /)
+ REAL(8), PARAMETER:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /)
TYPE, PUBLIC, EXTENDS(meshNode):: meshNodeCyl
!Element coordinates
@@ -284,7 +284,7 @@ MODULE moduleMeshCyl
END SUBROUTINE areaQuad
!Computes element functions in point xi
- PURE FUNCTION fPsiQuad(xi) RESULT(fpsi)
+ PURE FUNCTION fPsiQuad(xi) RESULT(fPsi)
IMPLICIT NONE
REAL(8),INTENT(in):: xi(1:3)
@@ -371,6 +371,7 @@ MODULE moduleMeshCyl
INTEGER:: l, m
ke=0.D0
+ xi=0.D0
!Start 2D Gauss Quad Integral
DO l=1, 3
xi(2) = corQuad(l)
@@ -489,7 +490,8 @@ MODULE moduleMeshCyl
CLASS(meshVolCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: dPsi(1:2,1:4)
- REAL(8):: j(1:2,1:2), detJ
+ REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates
+ REAL(8):: invJ(1:2,1:2), detJ
REAL(8):: phi(1:4)
REAL(8):: EF(1:3)
@@ -498,11 +500,12 @@ MODULE moduleMeshCyl
self%n3%emData%phi, &
self%n4%emData%phi /)
- dPsi = self%dPsi(xi)
- detJ = self%detJac(xi,dPsi)
- j = self%invJac(xi,dPsi)
- EF(1) = -DOT_PRODUCT(dPsi(1,:), phi)*j(2,2)/detJ
- EF(2) = -DOT_PRODUCT(dPsi(2,:), phi)*j(1,1)/detJ
+ dPsi = self%dPsi(xi)
+ detJ = self%detJac(xi,dPsi)
+ invJ = self%invJac(xi,dPsi)
+ dPsiR = MATMUL(invJ, dPsi)/detJ
+ EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi)
+ EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi)
EF(3) = 0.D0
END FUNCTION gatherEFQuad
@@ -526,7 +529,7 @@ MODULE moduleMeshCyl
CLASS(meshVolCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
- REAL(8):: xO(1:3), detJ, j(1:2,1:2), f(1:2)
+ REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2)
REAL(8):: dPsi(1:2,1:4), fPsi(1:4)
REAL(8):: conv
@@ -536,12 +539,12 @@ MODULE moduleMeshCyl
DO WHILE(conv>1.D-4)
dPsi = self%dPsi(xO)
- j = self%invJac(xO, dPsi)
+ invJ = self%invJac(xO, dPsi)
fPsi = self%fPsi(xO)
f(1) = DOT_PRODUCT(fPsi,self%z)-r(1)
f(2) = DOT_PRODUCT(fPsi,self%r)-r(2)
detJ = self%detJac(xO,dPsi)
- xN(1:2)=xO(1:2) - MATMUL(j, f)/detJ
+ xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ
conv=MAXVAL(DABS(xN-xO),1)
xO=xN
@@ -615,12 +618,39 @@ MODULE moduleMeshCyl
CLASS(meshVolCylTria), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
+ REAL(8):: ldCorner(1:3) !left down (ld) corner of square enclousing the triangle
+ INTEGER:: minP
+ INTEGER:: pOrder(1:3)
REAL(8), DIMENSION(1:3):: r1, r2, r3
+ !Assign node index
self%n = n
- self%n1 => mesh%nodes(p(1))%obj
- self%n2 => mesh%nodes(p(2))%obj
- self%n3 => mesh%nodes(p(3))%obj
+
+ !Find the first node
+ r1 = mesh%nodes(p(1))%obj%getCoordinates()
+ r2 = mesh%nodes(p(2))%obj%getCoordinates()
+ r3 = mesh%nodes(p(3))%obj%getCoordinates()
+ ldCorner = (/ MINVAL((/r1(1), r2(1), r3(1)/), 1), &
+ MINVAL((/r1(2), r2(2), r3(2)/), 1), &
+ 0.D0 /)
+ minP = MINLOC((/ NORM2(r1-ldCorner), NORM2(r2-ldCorner), NORM2(r3-ldCorner) /), 1)
+ pOrder(1) = p(minP)
+ SELECT CASE(minP)
+ CASE(1)
+ pOrder(2) = p(2)
+ pOrder(3) = p(3)
+ CASE(2)
+ pOrder(2) = p(3)
+ pOrder(3) = p(1)
+ CASE(3)
+ pOrder(2) = p(1)
+ pOrder(3) = p(2)
+ END SELECT
+
+ !Assign nodes to element
+ self%n1 => mesh%nodes(pOrder(1))%obj
+ self%n2 => mesh%nodes(pOrder(2))%obj
+ self%n3 => mesh%nodes(pOrder(3))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
@@ -653,7 +683,7 @@ MODULE moduleMeshCyl
self%arNodes = 0.D0
!2D 1 point Gauss Quad Integral
xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
- detJ = self%detJac(xi)*8.D0*PI
+ detJ = self%detJac(xi)*PI !2PI*1/2
fPsi = self%fPsi(xi)
r = DOT_PRODUCT(fPsi,self%r)
self%volume = r*detJ
@@ -697,7 +727,7 @@ MODULE moduleMeshCyl
REAL(8), INTENT(in):: xi2
REAL(8):: dPsiXi1(1:3)
- dPsiXi1(1) = -xi2
+ dPsiXi1(1) = -1.D0
dPsiXi1(2) = 1.D0
dPsiXi1(3) = 0.D0
@@ -710,7 +740,7 @@ MODULE moduleMeshCyl
REAL(8), INTENT(in):: xi1
REAL(8):: dPsiXi2(1:3)
- dPsiXi2(1) = -xi1
+ dPsiXi2(1) = -1.D0
dPsiXi2(2) = 0.D0
dPsiXi2(3) = 1.D0
@@ -739,20 +769,21 @@ MODULE moduleMeshCyl
REAL(8):: r, xi(1:3)
REAL(8):: fPsi(1:3), dPsi(1:2,1:3)
REAL(8):: ke(1:3,1:3)
- REAL(8):: j(1:2,1:2), detJ
+ REAL(8):: invJ(1:2,1:2), detJ
INTEGER:: l
ke=0.D0
+ xi=0.D0
!Start 2D Gauss Quad Integral
- DO l=1, 3
+ DO l=1, 4
xi(1) = xi1Tria(l)
xi(2) = xi2Tria(l)
dPsi = self%dPsi(xi)
detJ = self%detJac(xi,dPsi)
- j = self%invJac(xi,dPsi)
+ invJ = self%invJac(xi,dPsi)
fPsi = self%fPsi(xi)
r = DOT_PRODUCT(fPsi,self%r)
- ke = ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*wTria(l)/detJ
+ ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
END DO
ke = ke*2.D0*PI
@@ -776,7 +807,7 @@ MODULE moduleMeshCyl
localF = 0.D0
xi = 0.D0
!Start 2D Gauss Quad Integral
- DO l=1, 3
+ DO l=1, 4
xi(1) = xi1Tria(l)
xi(2) = xi2Tria(l)
detJ = self%detJac(xi)
@@ -852,7 +883,8 @@ MODULE moduleMeshCyl
CLASS(meshVolCylTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: dPsi(1:2,1:3)
- REAL(8):: j(1:2,1:2), detJ
+ REAL(8):: dPsiR(1:2,1:3)
+ REAL(8):: invJ(1:2,1:2), detJ
REAL(8):: phi(1:3)
REAL(8):: EF(1:3)
@@ -860,16 +892,17 @@ MODULE moduleMeshCyl
self%n2%emData%phi, &
self%n3%emData%phi/)
- dPsi = self%dPsi(xi)
- detJ = self%detJac(xi,dPsi)
- j = self%invJac(xi,dPsi)
- EF(1) = -DOT_PRODUCT(dPsi(1,:), phi)*j(2,2)/detJ
- EF(2) = -DOT_PRODUCT(dPsi(2,:), phi)*j(1,1)/detJ
+ dPsi = self%dPsi(xi)
+ detJ = self%detJac(xi,dPsi)
+ invJ = self%invJac(xi,dPsi)
+ dPsiR = MATMUL(invJ, dPsi)/detJ
+ EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi)
+ EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi)
EF(3) = 0.D0
END FUNCTION gatherEFTria
- !Gets nodes from triangular element
+ !Gets node indexes from triangular element
PURE FUNCTION getNodesTria(self) RESULT(n)
IMPLICIT NONE
@@ -882,33 +915,23 @@ MODULE moduleMeshCyl
END FUNCTION getNodesTria
!Transforms physical coordinates to element coordinates
- PURE FUNCTION phy2logTria(self,r) RESULT(xN)
+ PURE FUNCTION phy2logTria(self,r) RESULT(xi)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
- REAL(8):: xN(1:3)
- REAL(8):: xO(1:3), detJ, j(1:2,1:2), f(1:2)
- REAL(8):: dPsi(1:2,1:3), fPsi(1:3)
- REAL(8):: conv
+ REAL(8):: xi(1:3)
+ REAL(8):: invJ(1:2,1:2), detJ
+ REAL(8):: deltaR(1:2)
+ REAL(8):: dPsi(1:2,1:3)
- !Iterative newton method to transform coordinates
- !TODO: Try to find a direct alternative
- conv=1.D0
- xO=(/1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
-
- DO WHILE(conv>1.D-4)
- dPsi = self%dPsi(xO)
- j = self%invJac(xO, dPsi)
- fPsi = self%fPsi(xO)
- f(1) = DOT_PRODUCT(fPsi,self%z)-r(1)
- f(2) = DOT_PRODUCT(fPsi,self%r)-r(2)
- detJ = self%detJac(xO,dPsi)
- xN(1:2)=xO(1:2) - MATMUL(j, f)/detJ
- conv=MAXVAL(DABS(xN-xO),1)
- xO=xN
-
- END DO
+ !Direct method to convert coordinates
+ xi = (/ 0.D0, 0.D0, 0.D0 /) !Irrelevant, required for input
+ deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /)
+ dPsi = self%dPsi(xi)
+ invJ = self%invJac(xi, dPsi)
+ detJ = self%detJac(xi, dPsi)
+ xi(1:2) = MATMUL(invJ,deltaR)/detJ
END FUNCTION phy2logTria
@@ -969,7 +992,7 @@ MODULE moduleMeshCyl
CLASS(meshVolCyl), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
- REAL(8), INTENT(in), OPTIONAL:: dPsi_in(:,:)
+ REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dJ
REAL(8):: dz(1:2), dr(1:2)
@@ -993,8 +1016,8 @@ MODULE moduleMeshCyl
CLASS(meshVolCyl), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
- REAL(8), INTENT(in), OPTIONAL:: dPsi_in(:,:)
- REAL(8):: dPsi(1:2,1:4)
+ REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
+ REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dz(1:2), dr(1:2)
REAL(8):: invJ(1:2,1:2)
@@ -1033,7 +1056,7 @@ MODULE moduleMeshCyl
!
! END IF
part%e_p = self%n
- part%xi = xi
+ part%xi = xi
part%n_in = .TRUE.
!Assign particle to listPart_in
CALL OMP_SET_LOCK(self%lock)
diff --git a/src/modules/moduleMeshCylRead.f95 b/src/modules/moduleMeshCylRead.f95
index 9a0f5ff..3965afb 100644
--- a/src/modules/moduleMeshCylRead.f95
+++ b/src/modules/moduleMeshCylRead.f95
@@ -78,7 +78,6 @@ MODULE moduleMeshCylRead
BACKSPACE(10)
END DO
- !TODO: symplify to read all elements independently if they are edges or vols
!Reads edges
DO e=1, self%numEdges
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
@@ -97,16 +96,30 @@ MODULE moduleMeshCylRead
END SELECT
- CALL self%edges(e)%obj%init(n, p, bt, boundaryType)
+ CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType)
END DO
!Read and initialize volumes
- !TODO: Extend to triangular elements
DO e=1, self%numVols
- READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
- ALLOCATE(meshVolCylQuad:: self%vols(e)%obj)
- CALL self%vols(e)%obj%init(n - self%numEdges, p)
+ READ(10,*) n, elemType
+ BACKSPACE(10)
+
+ SELECT CASE(elemType)
+ CASE (2)
+ !Triangular element
+ READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
+ ALLOCATE(meshVolCylTria:: self%vols(e)%obj)
+ CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3))
+
+ CASE (3)
+ !Quadrilateral element
+ READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
+ ALLOCATE(meshVolCylQuad:: self%vols(e)%obj)
+ CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
+
+ END SELECT
+
END DO
@@ -115,10 +128,12 @@ MODULE moduleMeshCylRead
!Build connectivity between elements
DO e = 1, self%numVols
!Connectivity between volumes
- DO et = 1, self%numVols
- CALL connected(self%vols(e)%obj, self%vols(et)%obj)
+ IF (e /= et) THEN
+ DO et = 1, self%numVols
+ CALL connected(self%vols(e)%obj, self%vols(et)%obj)
- END DO
+ END DO
+ END IF
!Connectivity between vols and edges
DO et = 1, self%numEdges
@@ -133,19 +148,38 @@ MODULE moduleMeshCylRead
END SUBROUTINE
+ !Selects type of elements to build connection
SUBROUTINE connectedVolVol(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB
-
SELECT TYPE(elemA)
- TYPE IS (meshVolCylQuad)
+ TYPE IS(meshVolCylQuad)
+ !Element A is a quadrilateral
SELECT TYPE(elemB)
TYPE IS(meshVolCylQuad)
+ !Element B is a quadrilateral
CALL connectedQuadQuad(elemA, elemB)
+ TYPE IS(meshVolCylTria)
+ !Element B is a triangle
+ CALL connectedQuadTria(elemA, elemB)
+
+ END SELECT
+
+ TYPE IS(meshVolCylTria)
+ !Element A is a Triangle
+ SELECT TYPE(elemB)
+ TYPE IS(meshVolCylQuad)
+ !Element B is a quadrilateral
+ CALL connectedQuadTria(elemB, elemA)
+
+ TYPE IS(meshVolCylTria)
+ !Element B is a triangle
+ CALL connectedTriaTria(elemA, elemB)
+
END SELECT
END SELECT
@@ -158,12 +192,17 @@ MODULE moduleMeshCylRead
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
- SELECT TYPE(elemA)
- TYPE IS (meshVolCylQuad)
- SELECT TYPE(elemB)
- CLASS IS(meshEdgeCyl)
+ SELECT TYPE(elemB)
+ CLASS IS(meshEdgeCyl)
+ SELECT TYPE(elemA)
+ TYPE IS(meshVolCylQuad)
+ !Element A is a quadrilateral
CALL connectedQuadEdge(elemA, elemB)
+ TYPE IS(meshVolCylTria)
+ !Element A is a triangle
+ CALL connectedTriaEdge(elemA, elemB)
+
END SELECT
END SELECT
@@ -176,44 +215,208 @@ MODULE moduleMeshCylRead
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemB
- !Check direction 1
- IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
- elemA%n1%n == elemB%n4%n .AND. &
- elemA%n2%n == elemB%n3%n) THEN
+ !Check direction 1
+ IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+ elemA%n1%n == elemB%n4%n .AND. &
+ elemA%n2%n == elemB%n3%n) THEN
+ elemA%e1 => elemB
+ elemB%e3 => elemA
+
+ END IF
+
+ !Check direction 2
+ IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+ elemA%n2%n == elemB%n1%n .AND. &
+ elemA%n3%n == elemB%n4%n) THEN
+ elemA%e2 => elemB
+ elemB%e4 => elemA
+
+ END IF
+
+ !Check direction 3
+ IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
+ elemA%n3%n == elemB%n2%n .AND. &
+ elemA%n4%n == elemB%n1%n) THEN
+ elemA%e3 => elemB
+ elemB%e1 => elemA
+
+ END IF
+
+ !Check direction 4
+ IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
+ elemA%n4%n == elemB%n3%n .AND. &
+ elemA%n1%n == elemB%n2%n) THEN
+ elemA%e4 => elemB
+ elemB%e2 => elemA
+
+ END IF
+
+ END SUBROUTINE connectedQuadQuad
+
+ SUBROUTINE connectedQuadTria(elemA, elemB)
+ IMPLICIT NONE
+
+ CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
+ CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB
+
+ !Check direction 1
+ IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+ IF (elemA%n1%n == elemB%n1%n .AND. &
+ elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
+ ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
+ elemA%n2%n == elemB%n2%n) THEN
+ elemA%e1 => elemB
+ elemB%e2 => elemA
+
+ ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+ elemA%n2%n == elemB%n1%n) THEN
+ elemA%e1 => elemB
+ elemB%e1 => elemA
+
END IF
- !Check direction 2
- IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
- elemA%n2%n == elemB%n1%n .AND. &
- elemA%n3%n == elemB%n4%n) THEN
+ END IF
+
+ !Check direction 2
+ IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+ IF (elemA%n2%n == elemB%n1%n .AND. &
+ elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
- elemB%e4 => elemA
+ elemB%e3 => elemA
+
+ ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
+ elemA%n3%n == elemB%n2%n) THEN
+ elemA%e2 => elemB
+ elemB%e2 => elemA
+
+ ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+ elemA%n3%n == elemB%n1%n) THEN
+ elemA%e2 => elemB
+ elemB%e1 => elemA
END IF
- !Check direction 3
- IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
- elemA%n3%n == elemB%n2%n .AND. &
- elemA%n4%n == elemB%n1%n) THEN
+ END IF
+
+ !Check direction 3
+ IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+ IF (elemA%n3%n == elemB%n1%n .AND. &
+ elemA%n4%n == elemB%n3%n) THEN
+ elemA%e3 => elemB
+ elemB%e3 => elemA
+
+ ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
+ elemA%n4%n == elemB%n2%n) THEN
+ elemA%e3 => elemB
+ elemB%e2 => elemA
+
+ ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+ elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
+ END IF
+
+ END IF
+
+ !Check direction 4
+ IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+ IF (elemA%n4%n == elemB%n1%n .AND. &
+ elemA%n1%n == elemB%n3%n) THEN
+ elemA%e4 => elemB
+ elemB%e3 => elemA
+
+ ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
+ elemA%n1%n == elemB%n2%n) THEN
+ elemA%e4 => elemB
+ elemB%e2 => elemA
+
+ ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
+ elemA%n1%n == elemB%n1%n) THEN
+ elemA%e4 => elemB
+ elemB%e1 => elemA
END IF
- !Check direction 4
- IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
- elemA%n4%n == elemB%n3%n .AND. &
- elemA%n1%n == elemB%n2%n) THEN
- elemA%e4 => elemB
+ END IF
+
+ END SUBROUTINE connectedQuadTria
+
+ SUBROUTINE connectedTriaTria(elemA, elemB)
+ IMPLICIT NONE
+
+ CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
+ CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB
+
+ !Check direction 1
+ IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+ IF (elemA%n1%n == elemB%n1%n .AND. &
+ elemA%n2%n == elemB%n3%n) THEN
+ elemA%e1 => elemB
+ elemB%e3 => elemA
+
+ ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+ elemA%n2%n == elemB%n1%n) THEN
+ elemA%e1 => elemB
+ elemB%e1 => elemA
+
+ ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
+ elemA%n2%n == elemB%n2%n) THEN
+ elemA%e1 => elemB
elemB%e2 => elemA
END IF
- END SUBROUTINE connectedQuadQuad
+ END IF
+
+ !Check direction 2
+ IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+ IF (elemA%n2%n == elemB%n1%n .AND. &
+ elemA%n3%n == elemB%n3%n) THEN
+ elemA%e2 => elemB
+ elemB%e3 => elemA
+
+ ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+ elemA%n3%n == elemB%n1%n) THEN
+ elemA%e2 => elemB
+ elemB%e1 => elemA
+
+ ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
+ elemA%n3%n == elemB%n2%n) THEN
+ elemA%e2 => elemB
+ elemB%e2 => elemA
+
+ END IF
+
+
+ END IF
+
+ !Check direction 3
+ IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+ IF (elemA%n3%n == elemB%n1%n .AND. &
+ elemA%n1%n == elemB%n3%n) THEN
+ elemA%e3 => elemB
+ elemB%e3 => elemA
+
+ ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+ elemA%n1%n == elemB%n1%n) THEN
+ elemA%e3 => elemB
+ elemB%e1 => elemA
+
+ ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
+ elemA%n1%n == elemB%n2%n) THEN
+ elemA%e3 => elemB
+ elemB%e2 => elemA
+
+ END IF
+
+
+ END IF
+
+ END SUBROUTINE connectedTriaTria
SUBROUTINE connectedQuadEdge(elemA, elemB)
IMPLICIT NONE
@@ -287,6 +490,62 @@ MODULE moduleMeshCylRead
END SUBROUTINE connectedQuadEdge
+ SUBROUTINE connectedTriaEdge(elemA, elemB)
+ IMPLICIT NONE
+
+ CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
+ CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB
+
+ !Check direction 1
+ IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+ IF (elemA%n1%n == elemB%n1%n .AND. &
+ elemA%n2%n == elemB%n2%n) THEN
+ elemA%e1 => elemB
+ elemB%e2 => elemA
+
+ ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+ elemA%n2%n == elemB%n1%n) THEN
+ elemA%e1 => elemB
+ elemB%e1 => elemA
+
+ END IF
+
+ END IF
+
+ !Check direction 2
+ IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+ IF (elemA%n2%n == elemB%n1%n .AND. &
+ elemA%n3%n == elemB%n2%n) THEN
+ elemA%e2 => elemB
+ elemB%e2 => elemA
+
+ ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+ elemA%n3%n == elemB%n1%n) THEN
+ elemA%e2 => elemB
+ elemB%e1 => elemA
+
+ END IF
+
+ END IF
+
+ !Check direction 3
+ IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+ IF (elemA%n3%n == elemB%n1%n .AND. &
+ elemA%n1%n == elemB%n2%n) THEN
+ elemA%e3 => elemB
+ elemB%e2 => elemA
+
+ ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+ elemA%n1%n == elemB%n1%n) THEN
+ elemA%e3 => elemB
+ elemB%e1 => elemA
+
+ END IF
+
+ END IF
+
+ END SUBROUTINE connectedTriaEdge
+
SUBROUTINE constructGlobalK(K, elem)
IMPLICIT NONE
@@ -305,6 +564,13 @@ MODULE moduleMeshCylRead
n = (/ elem%n1%n, elem%n2%n, &
elem%n3%n, elem%n4%n /)
+ TYPE IS(meshVolCylTria)
+ nNodes = 3
+ ALLOCATE(localK(1:nNodes,1:nNodes))
+ localK = elem%elemK()
+ ALLOCATE(n(1:nNodes))
+ n = (/ elem%n1%n, elem%n2%n, elem%n3%n /)
+
CLASS DEFAULT
nNodes = 0
ALLOCATE(localK(1:1, 1:1))
@@ -461,6 +727,7 @@ MODULE moduleMeshCylRead
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
+ REAL(8):: xi(1:3)
IF (emOutput) THEN
time = DBLE(t)*tau*ti_ref
@@ -496,7 +763,13 @@ MODULE moduleMeshCylRead
WRITE(20, *) 3
WRITE(20, *) self%numVols
DO e=1, self%numVols
- WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF((/0.D0, 0.D0, 0.D0/))*EF_ref
+ SELECT TYPE(elem=>self%vols(e)%obj)
+ TYPE IS(meshVolCylQuad)
+ xi = (/ 0.D0, 0.D0, 0.D0 /)
+ TYPE IS(meshVolCylTria)
+ xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
+ END SELECT
+ WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
END DO
WRITE(20, "(A)") '$EndElementData'
CLOSE(20)