Denis Geraads

Nikolaï Paris

Sofia Spassov

Risto Garevski

Skopje

New specimens of Propotamochoerus (Suidae, Mammalia) from the late Miocene of the Balkans

Keywords: late Miocene, Turolian, Aegean region, Macedonia, Bulgaria, Mammalia, Suidae, Propotamochoerus 1. Introduction indet., Deinotherium gigantissimum, Ancylotherium pentelicum, Acerorhinus sp., Hipparion gr Bohlinia ? sp., Gazella sp., Prostrepsiceros sp., cf. "Ovis" kuhlmanni, cf, Tragoportax sp. 2. Systematic study Propotamochoerus hysudricus (STEHLIN, 1900)

We describe new remains of the small Suidae Propotamochoerus from the late Miocene of Macedonia and Bulgaria, which greatly increase the European sample. They are clearly distinct from both the Vallesian P. palaeochoerus and the Pliocene P. provincialis, but cannot be satisfactorily referred to the Asian species P. hysudricus or P. hyotherioides. It is likely that these Balkan specimens belong to a species of their own, probably distinct from the P. palaeochoerus -P. provincialis lineage.

it is clear that the talonid consisted of two tubercles, the lingual one (hexaconid) being only slightly smaller than the main one (pentaconid).

The material from Kalimantsi consists of an isolated P2 and a M3. On P2 (Fig. 2E) the protocone is quite strong, giving the tooth a high W/L ratio, but it is too narrow for a P3. The M3 (Fig. 2D) is slightly narrower anteriorly than the Macedonian specimens, but otherwise similar, with the talon in line with the lingual cusps.

Comparisons

The suids from Macedonia and Bulgaria can be compared to several taxa from the late Miocene and early Pliocene, which are smaller than Microstonyx, the common European suid of this period, but larger than Sus, and are usually assigned to the genus Propotamochoerus.

Propotamochoerus palaeochoerus (KAUP, 1833) is best known from the Vallesian of Central Europe (MOTTL 1966;HÜNERMANN 1968;SCHMIDT-KITTLER 1971) and Ukraine (VAN DER MADE et al. 1999). A record of earlier age is Montréjeau, France (HÜNERMANN 1968), but reports from younger localities are doubtful. A suid from Maramena, Serres, N. Greece, was described under this name by HELLMUND (1995), and a maxilla from Salihpaşalar, Turkey, is also labelled as such in the Museum of the General Directorate of Mineral Research and Exploration, Ankara (MTA-2383). The material from Alcoy, Spain, previously assigned to this species has recently been re-assigned to Sus arvernensis (MONTOYA et al. 2006). Propotamochoerus palaeochoerus was also listed from Tataros in Romania by RABEDER (1985); the material is too poor for specific identification but the genus is probably correct.

Propotamochoerus hysudricus (STEHLIN, 1900) from the Indian subcontinent has a complex synonymy list (reviews in PICKFORD 1988;BONIS & BOUVRAIN 1996), but is unfortunately poorly illustrated (PILGRIM 1926;PICKFORD 1988). A mandible of a similar form has been described by BONIS & BOUVRAIN (1996) from the early Turolian of Ravin des Zouaves n°5 (RZO), and another has been discovered by E. TSOUKALA at Thermopigi, both in Northern Greece.

Propotamochoerus hyotherioides (SCHLOSSER, 1903) is mostly known by isolated teeth from Lufeng, China (VAN DER MADE & HAN, 1994) but PEARSON (1928) also illustrated a palate and a mandible from Loc. 49. This Chinese species was included in Hippopotamodon by PICKFORD & LIU (2001), and there are indeed a number of similarities between it and both species of this genus, H. sivalense from the Middle Siwaliks of Pakistan, and H. antiquus from the early Late Miocene of Western Eurasia, which are significantly larger, however. Pending revision of the whole group, we follow common practice in referring hyotherioides to Propotamochoerus. THENIUS (1950) referred to the same species a mandible and associated maxilla from Samos that BONIS & BOUVRAIN (1996) compared instead to the RZO form. MITZOPOULOS (1947) also reported a maxilla fragment from Achladi in Euboea, while GALLAI & ROOK (2006) referred to this species or the preceding one a few specimens from Sicily.

Propotamochoerus provincialis (GERVAIS, 1852) has been reported from several latest Miocene and early Pliocene European sites, but is best known from the early Pliocene of Montpellier; the material has never been published in detail, but we were able to study a good sample in "Université Claude Bernard" and "Centre de Conservation et d' Etude des Collections", Lyon. The species is also present at Malusteni, Romania (SIMIONESCU 1930) but the specimens from Kvabebi, Georgia (VEKUA 1972;VEKUA et al. 2005) are clearly not of this species, as noted by AZZAROLI (1975).

P1: In the Late Miocene, this tooth is bi-rooted, bi-cuspidate, and low, as at Vozarci. On the male specimen of P. palaeochoerus described by MOTTL (1966), the diastema separating this tooth from P2 is much longer than in all other forms. Specimen FSL-40073 of P. provincialis from Montpellier displays only the root of a monoradiculate vestigial P1, almost contacting both the canine and P2; this tooth was therefore much more reduced than at Vozarci. P2: In P. palaeochoerus, and in P. hyotherioides from Lufeng, P2 is distinctly narrower anteriorly than posteriorly, and the protocone is better-marked than at Vozarci or at Samos, where the tooth is of more regular width. The tooth from Kalimantsi has a very strong protocone, and is the broadest of all Propotamochoerus P2s, but there is no other taxon to which it could be referred. At Maramena P2 is longer than P3, as at Vozarci, whereas they are about as long in other sites, but in P. provincialis, instead, this tooth is smaller than P3.

P3: The protocone is variable in P. palaeochoerus, but in this species as well as in the other ones, it is always larger than at Samos. The main cusp is usually double, except in specimens from the Indian subcontinent assigned to P. hysudricus, where the main cusp is stout and rounded, more like at Maramena and in P. provincialis. P4: This tooth has a rather constant morphology. The paracone is usually slightly larger than the metacone, but the difference in size between these cusps is more important, and the protocone slightly shifted mesially, in the specimen from Gravitelli, in P. hyotherioides from Lufeng, and especially in P. provincialis.

M3:

The size of the talon (pentacone of VAN DER MADE, 1996) is not very variable, but in P. palaeochoerus it is situated centrally, so that the tooth is almost symmetrical. One exception is the specimen from Grytsiv. In younger forms, instead, the pentacone is lingually shifted, and almost in line with the protocone and "hypocone" (tetracone of VAN DER MADE, 1996), and the tooth is no longer symmetrical, but more triangular. This might be connected with the overall trend towards development of the hexaconid on m3 (see below).

Lower canine: The mandible from Vozarci seems to have canines of scrofa-type, as in all other taxa. HÜNERMANN (1968) showed that both the size and detailed morphology are rather variable within P. palaeochoerus. p1: In the Indian P. hysudricus, this tooth is inserted immediately behind the canine. In the Vozarci mandible, it is inserted less than 1 cm behind it, but a longer diastema separates it from p2. The condition is very similar in the RZO and Samos mandibles, and in P. provincialis, but the diastema between p1 and p2 is shorter, if not absent, in P. hysudricus from India. p2: This tooth displays little variation, except that those from RZO are large. A premolar from Tataros, Romania, is definitely too small to be a p3, but is even smaller and broader than those of P. palaeochoerus; it is probably best assigned to this latter taxon. p3: In the Macedonian specimens, this tooth is broad (Fig. 3) but oval-shaped, and more similar to p2 than to p4; the same is true of p3s from RZO. In P. palaeochoerus, p3 is also ovalshaped, but quite narrow. The p3s of the Indian P. hysudricus are always short, in contrast to the P. cf. hysudricus from RZO. p4: The p4s of P. palaeochoerus (HLMD Din-3, Din-48, Din-94) are quite distinct from those of the Balkans by their strong protruding Innenhügel and by their talonid remaining long isolated from the trigonid. VAN DER MADE & HAN (1994) illustrated only unworn or little worn p4s of P. hyotherioides from Lufeng, but they look more similar to those of P. palaeochoerus than to those of the Balkan form, because the main cuspid is clearly subdivided, and the grooves limiting the talonid more deeply incised. Most of the p4s of P. provincialis that we have seen are not much worn but, even on the oldest specimen (FSL 40884), the talonid remains well separated from the trigonid, which has a rather strong metaconid. By contrast, on the p4 of the mandible from Ravin des Zouaves referred to P. cf. hysudricus, the dentine island is continuous as at Vozarci, although the tooth is not much worn (BONIS & BOUVRAIN 1996, Fig. 9).

If we compare the length proportions, the pattern is clearer than for p3 (Fig. 4). All specimens assigned to P. palaeochoerus (except one aberrant mandible from Hennersdorf: HELLMUND 1995) have long p4s (as noted by VAN DER MADE et al. 1999), whereas the p4/m3 ratio is smaller and more or less the same in all other forms, including those from Vozarci, RZO and Lufeng, that have long p3s.

The m1s and m2s, often much worn, do not display significant differential features. m3: An important issue is the significance of the variations in the composition of the distal complex. There is some variation at Vozarci, and it is well-known that Plio-Pleistocene African suids display a great range of variation in this structure, even within a single stratigraphic level. On the other hand, VAN DER MADE et al. insisted (VAN DER MADE et al. 1999: 286) on the occurrence of two pillars in the distal lobe of the RZO specimen, and suggested (VAN DER MADE et al. 1999: 290) that this has taxonomic significance, the hexaconid being absent in P. palaeochoerus and P. hysudricus. However, the sample of P. provincialis from Montpellier clearly demonstrates that intra-specific variation may be quite large. In some specimens from this locality (FSL 40155, FSL 40161), the distal complex consists almost only of the pentaconid, with only minor accessory tubercles; others have a stronger hexaconid (FSL 40070,FSL 40072); others (FSL 40175,FSL 40881) have an almost symmetrical distal complex, with the hexaconid even larger than the pentaconid (but no hint of a heptaconid). Thus, the variation at Montpellier is even larger than at Vozarci. It is therefore clear that the occurrence of one or two cuspid in the distal complex of m3 is insufficient to support a taxonomic distinction. What remains true, however, is that there is a clear trend, from the Vallesian to the Ruscinian, for the distal complex to become more often bicuspid.

We performed a principal component analysis on the most significant lower teeth, excluding p1 and p2 because they are too rare, and m1-m2 because their measurements are often reduced by wear. Results, shown in Fig. 5, mostly confirm those of the bivariate plots.

Propotamochoerus provincialis is distinct by its large size. The variation of the Indian P. hysudricus is probably too large to be accommodated within a single species.

Propotamochoerus palaeochoerus from the Vallesian has long and broad p4s, and is well-distinct in this regard from all other forms, including the one from Maramena. The third axis of the PCA (inertia: 7.68 %) describes mostly the width of the premolars, and separates clearly the Macedonian specimens (but not so clearly the RZO form) from the Lufeng P. hyotherioides.

Discussion and conclusions

It is likely that all the forms from Kalnitsa and Vozarci (Macedonia), Kalimantsi (Bulgaria), RZO, Samos, Achladi, and Thermopigi (Greece), all belong to the same species, which remains relatively rare and co-exists, at Vozarci, Kalimantsi, RZO, and Thermopigi, with the large Microstonyx. From Salihpaşalar (Turkey), PICKFORD & ERTÜRK (1979) reported only Microstonyx, but the maxilla mentioned above is similar in size and morphology to Propotamochoerus, and likely belongs to this genus as well, rather than to Hippopotamodon, represented in Turkey by a larger and earlier species. The material from Maramena is probably of the same taxon; it lacks the features of P. palaeochoerus, especially the large p4 and symmetrical M3, but in spite of the relatively recent age of the locality, the third molars are smaller than in P. provincialis.

Propotamochoerus palaeochoerus is well-distinct from this Turolian form, and it may even be doubted that they belong to the same genus; thus, "Korynochoerus" palaeochoerus might indeed deserve a distinct generic name. The skull from Johnsdorf described by MOTTL (1966) has very large supra-canine flanges, and the premolars are quite distinct: P1 is small and separated from P2 by a long diastema, P2 is also short but much broader distally than mesially, P3 has a very large protocone giving the tooth a triangular outline, and P4 is shorter lingually than labially. On the whole, the premolars are significantly more reduced than at Vozarci. Furthermore, in P. palaeochoerus, the talon of M3 is centrally placed, instead of being almost in line with the protocone and hypocone. Besides these dental differences, compared with that of P. hysudricus, the skull and especially the snout are longer, the supra-canine flange is larger, the zygomatic arches more expanded anteriorly, and the temporal lines meet in the sagittal plane.

Propotamochoerus provincialis has slightly longer M3s and m3s than the Balkan form, P1 and P2 are reduced, the protocone of P4 is shifted mesially, p3 is broad distally, the talonid of p4 is well-distinct from the main cuspid, and the third lobe of m3 trends towards bicuspidy.

Interestingly, some features of this species are shared by P. palaeochoerus, but not by the Eastern Mediterranean form, and one could imagine that P. provincialis is an independent offshoot from the Vallesian form.

The comparison of the Balkan form with P. hysudricus and P. hyotherioides is less straightforward, because:

-both named species are based upon inadequate material and remain poorly defined; -as we have seen above, and as already noted by BONIS & BOUVRAIN (1996), the Indian material referred P. hysudricus probably includes more than one species; -cranial material is known only for "P. hysudricus".

Assignment of the Balkan form to either P. hysudricus or P. hyotherioides would not be very satisfactory, and it is likely that it belongs in fact to an as yet unnamed species, restricted to the Aegean region, but this remains to be further supported by cranial material. Its age remains at the present time poorly constrained. Ravin des Zouaves has been assigned to "zone" MN 11; the tooth from Kalimantsi is likely from the first half of the middle Turolian, but the faunas of the Republic of Macedonia have not yet been studied in detail, although an MN 11 or early MN 12 age is likely for Vozarci. We may assume, as a working hypothesis, that most Propotamochoerus localities of the Balkans and Turkey, except Maramena, belong to the first half of the Turolian.

The last species of the genus, P. provincialis, became extinct in the Middle Pliocene; it is at that time that a related genus, Kolpochoerus, is first recorded Africa, first by K. deheinzelini BRUNET & WHITE, 2001, a species which is still poorly known and of uncertain status (GERAADS, 2004), and later, at about 3.5 -3 Ma., by K. afarensis COOKE, 1978. A close relationship between both genera is generally accepted, and one may conceive that some species of Propotamochoerus gave rise to the African forms, which are more derived by their less straight cranial profile, backward retraction of the zygoma and longer muzzle, larger upper canine, and more complex P3 (with a mesio-lingual cusp which is always absent in Propotamochoerus). However, P. provincialis is an unlikely ancestor for Kolpochoerus because of its thick premolars, mesially shifted protocone on P4, and talonid of m3 more complex than in early Kolpochoerus, so that if the latter arose from a Eurasian immigrant, this must have happened in the latest Miocene, which would confirm K. deheinzelini as a member of this genus.

van Wetenschappen, 97: 27-82. VAN (HÜNERMANN 1968;HELLMUND 1995), the mean values of P. palaeochoerus from Maramena (HELLMUND 1995), the mean values of P. palaeochoerus from Rudabánya (FORTELIUS et al. 2005), two well-preserved specimens of P. palaeochoerus from Grytsiv (VAN DER MADE et al. 1999), the mandible of P. cf. hysudricus from Ravin des Zouaves 5 (BONIS & BOUVRAIN 1996), the mean values of P. hyotherioides from

  Figure 1: map of the Eastern Mediterranean, with the late Miocene Propotamochoerus localities.

Figure 2 :

 2 Figure 2: Propotamochoerus sp. Upper teeth are shown as if from the right side, lower teeth as if from the left side. A: upper right tooth row, Vozarci-271; B: upper left (reversed) M2-M3, Vozarci-895; C: upper right C-P3, Vozarci-565; D: upper right M3 from Kalimantsi; E: upper right P2 from Kalimantsi; F: lower left p3-m3, Vozarci-564; G: lower left p2-m3 from Kalnitsa; H: unnumbered mandible with left c-p1 and right p2-m3 (reversed), from Vozarci. Scale = 10 cm.

Fig. 3 :

 3 Fig. 3: Plot of p4 width vs. p3 width in Propotamochoerus. Mean values are shown by larger symbols.

Fig. 5 :

 5 Fig. 5: Plane 1-2 of the PCA on lower teeth measurements, using : well preserved tooth-rows and the mean values of P. hysudricus from India (PICKFORD 1988), the mean values of P. palaeochoerus from Wissberg (HÜNERMANN 1968), the mean values of a number of other Vallesian specimens of P. palaeochoerus from Central Europe (HÜNERMANN 1968; HELLMUND

5. Acknowledgements D. BERTHET (Centre de Conservation et d'Etude des Collections, Lyon), G. GRÜBER (Hessisches Landesmuseum, Darmstadt), M. IVANOV (Department of Geology, University of Sofia), and A. PRIEUR (Université Claude Bernard, Lyon), kindly gave us access to collections in their care. V. CODREA and M. GASPARIK provided data about the Tataros suid. Special thanks to L. HRISTOVA for her preliminary identifications of hipparions, to D. NAGEL for helpful comments, and to E. TSOUKALA for showing us her Thermopigi fossils. Remarks and suggestions by M. PICKFORD and an anonymous reviewer greatly improved the manuscript. All travel funds were provided by the CNRS, UPR 2147. 6. References AZZAROLI, A. (1975): Remarks on the Pliocene Suidae of Europe. -Zeitschrift für Säugetierkunde, 40: 355-367. BONIS, L. de & BOUVRAIN, G. (1996): Suidae du Miocène supérieur de Grèce. -Bulletin du Muséum national d'Histoire naturelle, 18: 107-132. BOUVRAIN, G. & BONIS, L. de (2007): Ruminants (Mammalia, Artiodactyla: Tragulidae, Cervidae, Bovidae) des gisements du Miocène supérieur (Turolien) de Dytiko (Grèce). -Annales de Paléontologie, 93: 121-147. BRUNET, M. & WHITE, T.D. (2001): Deux nouvelles espèces de Suini (Mammalia, Suidae) du continent africain (Ethiopie; Tchad). -Comptes-Rendus de l'Académie des Sciences, Sciences de la Terre et des Planètes, 332: 51-57. COOKE, H.B.S. (1978): Pliocene-Pleistocene Suidae from Hadar, Ethiopia. -Kirtlandia, 29: 1-63. FORSTEN, A. & GAREVSKI, R. (1989): Hipparions (Mammalia, Perissodactyla) from Macedonia, Yugoslavia. -Geologica macedonica, 3: 159-206. FORTELIUS, M., ARMOUR-CHELU, M., BERNOR, R. L. & FESSAHA, N. (2005): Systematics and Palaeobiology of the Rudabánya Suidae. -Palaeontographia Italica, 90: 259-278. GALLAI, G. & ROOK, L. (2006): Propotamochoerus sp. (Suidae, Mammalia) from the late Miocene of Gravitelli (Messina, Sicily, Italy) rediscovered. -Rivista Italiana di Paleontologia e Stratigrafia, 112: 317-321. GERAADS, D. (2004): New skulls of Kolpochoerus phacochoeroides (Suidae, Mammalia) from the late Pliocene of Ahl al Oughlam, Morocco. -Palaeontologia africana, 40: 69-83. GERAADS D., SPASSOV N. & KOVACHEV D. (2003): Palaeoreas lindermayeri (Wagner, 1848) from the upper Miocene of Bulgaria, and a revision of the species. -Géodiversitas, 25(2): 405-415. GERAADS, D., SPASSOV, N. & KOVACHEV D. (2005): Giraffidae (Artiodactyla, Mammalia) from the late Miocene of Kalimantsi and Hadjidimovo, South-Western Bulgaria. Geologica Balcanica, 35 (1-2): 11-18. GERVAIS, P. (1848-1852): Zoologie et Paléontologie françaises. -271 p.; Paris (A. Bertrand). HELLMUND, M. (1995): The Vertebrate Locality Maramena (Macedonia, Greece) at the Boundary Turolian-Ruscinian Boundary (Neogene). 13 Suidae (Artiodactyla, Mammalia). -Münchner Geowissenschaftliche Abhandlungen, 28: 143-156. HRISTOVA, L. & KOVACHEV, D. (2005): Notes on the Taxonomy of the Middle Turolian Hipparions from the Localities near Kalimantsi Village. Proceedings of the International Scientific Conference (8-11.06.2005). Faculty of Mathematics and Natural Sciences. South-Western University, Blagoevgrad, 180-184. (in Bulgarian, Engl. summary). HÜNERMANN, K. A. (1968): Die Suidae (Mammalia, Artiodactyla) aus den Dinotheriensanden (Unterpliozän= Pont) Rheinhessens (Südwestdeutschland). -Schweizerische Paläontologische Abhandlungen, 86: 1-96. KOSTOPOULOS, D. S., SPASSOV, N. & KOVACHEV, D. (2001): Contribution to the study of Microstonyx: evidence from Bulgaria and the SE European populations. -Geodiversitas, 23(3): 411-437. KOUFOS, G. D., SPASSOV, N. & KOVACHEV, D. (2003): Study of Mesopithecus (Primates, Cercopithecidae) from the late Miocene of Bulgaria. -Palaeontographica A, 269: 39-91. MARKOV, G. N. (2004). The fossil Proboscideans of Bulgaria. Unpubl. PhD Thesis, Sofia, 225 p. (in Bulgarian with English summary). MITZOPOULOS, M. K. (1947): Die Verbreitung der Pikermistufe auf der Insel Euböa. -Annales géologiques des Pays Helléniques, 1: 209-216. MONTOYA, P., GINSBURG, L., ALBERDI, M.-T., VAN DER MADE, J., MORALES, J. & SORIA, M.-D. (2006): Fossil large mammals from the early Pliocene locality of Alcoy (Spain) and their importance in biostratigraphy. -Geodiversitas 28(1): 137-173. MOTTL, M. (1966): Ein vollständiger Hyotherium palaeochoerus-Schädel aus dem Altpliozän (Pannon) Südost-Österreichs. -Mitteilungen des Museums für Bergbau, Geologie und Technik am Landesmuseum "Joanneum", 28: 73-103. NIKOLOV, I. (1985): Catalogue of the localities of tertiary Mammals in Bulgaria. -Paleontology, Stratigraphy and Lithology, 21: 43-62. PEARSON, H. S. (1928): Chinese fossil Suidae. -Palaeontologia Sinica, 5: 1-75. PICKFORD, M. (1988): Revision of the Miocene Suidae of the Indian Subcontinent. -Münchner Geowissenschaftliche Abhandlungen, 12: 1-95. PICKFORD, M. & ERTÜRK, Ç. (1979): Suidae and Tayassuidae from Turkey. -Bulletin of the Geological Society of Turkey, 22: 141-154. PICKFORD, M. & LIU LIPING (2001): Revision of the Miocene Suidae of Xiaolongtan (Kaiyuan), China. -Bolletino delle Società Paleontologica Italiana, 40: 275-283. PILGRIM, G. E. (1925): Presidential address to the geological section of the 12th Indian Science company. -Proceedings of the 12th Indian Science Congress: 200-218. [not seen; fide PICKFORD, 1988] PILGRIM, G. E. (1926): The fossil Suidae of India. -Memoirs of the Geological Survey of India. Palaeontologia Indica, 8(4): 1-65. RABEDER, G. (1985): Die Säugetiere des Pannonien. -In: PAPP, A., JAMBOR, A. & STEININGER, F. (Eds): Chronostratigraphie und Neostratotypen. Miozän der Zentralen Paratethys. M6 Pannonien, 440-463; Budapest (Ungarischer Akademie der Wissenschaften). SCHLOSSER, M. (1903): Die fossilen Säugetiere Chinas nebst einer Odontographie der rezenten Antilopen. -Abhandlungen der bayerischen Akademie der Wissenschaften, 7: 371-464. SCHMIDT-KITTLER, N. (1971): Die obermiozäne Fossillagerstätte Sandelzhausen, 3. Suidae (Artiodactyla, Mammalia). -Mitteilungen der bayerischen Staatssammlung für Paläontologie und historische Geologie, 11: 129-170. SIMIONESCU, I. (1930): Vertebratele Pliocene dela Malusteni (Covurlui). -Academia Romana-Publicatiunile Adamachi, 9: 83-148. SPASSOV, N., TZANKOV, TZ. & GERAADS, D. (2006): Late Neogene stratigraphy, biochronology, faunal diversity and environments of South-West Bulgaria (Struma River valley). -Géodiversitas, 28(3): 477-498. STEHLIN, H. (1900): Ueber die Geschichte des Suiden-Gebisses. -Mémoires Suisses de Paléontologie, 26: 1-527. THENIUS, E. (1950): Postpotamochoerus nov. subgen. hyotherioides aus dem Unterpliozän von Samos (Griechenland) und die Herkunft der Potamochoeren. -Sitzungsberichten der Österreichischen Akademie der Wissenschaften Mathematisch-Naturwissenschaftliche Klasse, 159: 25-36. VAN DER MADE, J. (1996): Listriodontinae (Suidae, Mammalia), their evolution, systematics and distribution in time and space. Contributions to Tertiary and Quaternary geology, 33: 3-254. VAN DER MADE, J. & HAN, D. (1994): Suoidea from the upper Miocene hominoid locality of Lufeng, Yunnan province, China. -Proceedings of the Koninklijke Nederlandse Akademie Lufeng (VAN DER MADE & HAN 1994) and the mandibles Vozarci-564 and Kalnitsa. Illustrative individuals are italicized. Mean values are represented by larger symbols. Addresses of the authors: DENIS GERAADS, UPR 2147 CNRS, 44 rue de l'Amiral Mouchez, 75014 Paris, France, dgeraads@ivry.cnrs.fr (corresponding author) NIKOLAÏ SPASSOV, National Museum of Natural History, Tsar Osvoboditel Blvd 1, 1000 Sofia, Bulgaria, nspassov@nmnh.bas.bg RISTO GAREVSKI, Macedonian Museum of Natural History, Ilinden Blvd 86, 91000 Skopje, Republic of Macedonia

Table 1 . Measurements of the upper teeth


Matania Ben-Artzi

Jean-Pierre Croisille
email: jean-pierre.croisille@univ-lorraine.fr

Dalia Fishelov
email: daliaf@afeka.ac.il

An embedded compact scheme for biharmonic problems in irregular domains

a Cartesian embedded finite difference scheme for biharmonic problems has been introduced. The design of the scheme relies on a 19-dimensional polynomial space. In this paper, we show how to simplify the implementation by introducing a directional decomposition of this space. The boundary is handled via a level-set approach. Numerical results for non convex domains demonstrate the fourth order accuracy of the scheme.

Introduction

Let Ω ⊆ R 2 be a convex domain. The problem considered here is the biharmonic problem subject to Dirichlet boundary conditions:

∆ 2 ψ(x) = f , x ∈ Ω , ψ = ∂ ψ ∂ n = 0, x ∈ ∂ Ω . ( 1
)
Our purpose is to calculate a high order accurate approximation to [START_REF] Abarbanel | Asymptotically Stable Fourth-Order Accurate Schemes for the Diffusion Equation on Complex Shapes[END_REF], by embedding Ω in a Cartesian grid. The main idea of the scheme was described in [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF]. Here we extend and elaborate on the presentation in [START_REF] Ben-Artzi | Navier-Stokes equations in planar domains[END_REF]Chap.11].

We consider the convex domain Ω as embedded in a large uniform grid of mesh size h. A grid point is a point Q i, j = (ih, jh) for i, j ∈ Z. Following common terminology, we use the term interior nodes for the grid points that lie inside Ω . We denote by Ω h the ensemble of these nodes, namely:

Ω h = Q i, j ∈ Ω , i, j ∈ Z . (2)
We split the set Ω h into two sets,

Ω h = Ω calc h ∪ Ω edge h
, as follows:

• Ω calc h = the set of calculated nodes. This set consists of those nodes that are located "well within" Ω , namely sufficiently far from the boundary ∂ Ω . In particular, if all diagonally neighboring nodes Q i±1, j±1 are in Ω h then Q i, j ∈ Ω calc h . Remark that by convexity all eight neighboring nodes are then in Ω h . However, it should be emphasized that even if not all its neighboring nodes are in Ω h , a node Q i, j can still be considered as "calculated" if it is not "too close" to the boundary, as we explain below. The approximate values at the calculated nodes are obtained by the proposed scheme.

• Ω edge h = the set of edge nodes. This set consists of those nodes (interior to Ω ) that are located "too close" to the boundary ∂ Ω . They differ from the calculated nodes in the sense that there are no approximate values associated with them. Their role is "geometric"; they serve in the determination of a set Ω bdry h of boundary nodes that are actually located on the boundary ∂ Ω , and carry the assigned boundary values.

• Observe that the set Ω bdry h consists of selected points on the boundary, and in general is not included in the underlying global grid Q i, j , i, j ∈ Z.

In Fig. 1 we designate the calculated nodes with black circles, whereas the edge nodes are designated by white circles.

The proposed scheme is a compact scheme, i.e. all approximate values of high order derivatives are related to values of a function ψ and its derivatives ψ x , ψ y at immediate neighbors. More specifically, given a node M 0 = Q i, j ∈ Ω h , we consider the eight surrounding nodes in the grid:

M 1 = Q i-1, j+1 , M 2 = Q i, j+1 , M 3 = Q i+1, j+1 , M 4 = Q i-1, j , M 5 = Q i+1, j , M 6 = Q i-1, j-1 , M 7 = Q i, j-1 , M 8 = Q i+1, j-1 .
If all the nine nodes M i are calculated nodes, namely, in Ω calc h , or coincide with a boundary point, which is part of the grid, we set M i = M i , i = 0, ..., 8, and continue with this regular stencil centered at M 0 . Otherwise, our goal is to replace the M i s that are not in Ω calc h by suitable M i s that are boundary points, namely, in Ω bdry h . The values of ψ, ψ x , ψ y at these points are all that is needed in order to calculate the various approximate derivatives at M 0 .

M 0 (ih, jh) Fig. 1 Embedding of an elliptical domain in a Cartesian grid. Calculated nodes are represented by black circles. Exterior points are represented by black squares. The points labelled with white circles represent edge points, i.e. interior points close to the boundary.

To describe this construction, suppose that M 0 ∈ Ω calc h is a calculated node, while for some 1 ≤ i ≤ 8, the neighboring node M i is either an edge node or an exterior node. Consider the calculated node designated by M 0 on Fig. 1. A zoom is shown on Fig. 2. The 8 points M i are the points on the square (4 corner points and 4 midedge points). Take the ray that emanates from M 0 and goes through M i . This ray must cross the boundary ∂ Ω at exactly one point since Ω is convex. We define the intersection point as M i .

The calculation of the approximate value to ∆ 2 ψ(M 0 ) relies on the data at M i rather than M i .

• The four neighbors M 1 , M 4 , M 6 and M 7 are other calculated nodes so we keep them, i.e. M i = M i , i = 1, 4, 6, 7. In particular, if we shift the coordinates of M 0 to (0, 0), we have for the coordinates of M i , i = 1, 4, 6, 7, the values h 1 = h 4 = h 6 = h 7 = h. • The other four neighbors M 2 , M 3 , M 5 and M 8 are either edge or exterior nodes so they are replaced by points on the boundary as described above.

We thus obtain M i , the actual points used in the calculation. Once the 8 points M i are determined and approximate values ψ,ψ x and ψ y are assigned to them, we can proceed to evaluate an approximate value for ∆ 2 ψ at the point M 0 . This is described in Section 2.

M 1 (-h 1 , h 1 ) M 2 (0, h 2 ) M 3 (h 3 , h 3 ) M 4 (-h 4 , 0) M 5 (h 5 , 0) M 6 (-h 6 , -h 6 ) M 7 (0, -h 7 ) M 8 (h 8 , -h 8 ) M 0 (0, 0) M 2 M 3 M 5 M 8

The discrete biharmonic ∆ h ψ operator

In this section we present our finite-difference scheme for the approximation of the biharmonic operator. Fig. 2 shows the stencil used for the approximation of ∆ 2 ψ at M 0 = (0, 0). The 8 points M k , 1 ≤ k ≤ 8 form an irregular stencil around M 0 . Each of the nine grid points M k carries three values: ψ, ψ x , ψ y . These are calculated values if

M k ∈ Ω calc h is a calculated node. If M k ∈ Ω bdry h
, then this point carries boundary data given by the boundary conditions. In order to approximate ∆ 2 ψ of a given smooth function ψ at M 0 we interpolate the data ψ, ψ x , ψ y on the stencil {M 0 , ..., M 8 } by a certain polynomial P M 0 of degree 6. The detailed construction of P M 0 (x, y) is carried out in Section 3. To handle the irregular stencil around M 0 we denote by h the vector of the step-sizes around M 0 , as in Figure 2:

h = [h 1 , ..., h 8 ] T . (3)
Once the polynomial P M 0 (x, y) is constructed, we replace the smooth function ψ by a discrete function ψ, defined only on the set of nodes

Ω calc h ∪ Ω bdry h . The dis- crete biharmonic operator ∆ 2 h ψ for the approximation of ∆ 2 ψ at M 0 = (0, 0) is then defined by ∆ 2 h ψ(M 0 ) = ∆ 2 P M 0 (0, 0), (4)
3 Calculating the interpolation polynomial P M 0 (x, y)

As mentioned above, our compact scheme for the biharmonic problem relies on an interpolation polynomial of degree six. Such a polynomial is constructed for every calculated point M i, j ∈ Ω calc h . This sixth-order polynomial is called P M 0 (x, y). It is of the form (where here and below the subscript M 0 is omitted) ,

P(x, y) = 19 ∑ i=1 a i l i (x, y), (5)
where the polynomials l i (x, y) are ((x, y) are shifted so that M 0 = (0, 0):

                           l 1 (x, y) = 1, l 2 (x, y) = x, l 3 (x, y) = x 2 , l 4 (x, y) = x 3 , l 5 (x, y) = x 4 , l 6 (x, y) = x 5 , l 7 (x, y) = y, l 8 (x, y) = y 2 , l 9 (x, y) = y 3 , l 10 (x, y) = y 4 , l 11 (x, y) = y 5 , l 12 (x, y) = xy, l 13 (x, y) = xy(x + y), l 14 (x, y) = xy(x -y), l 15 (x, y) = xy(x + y) 2 , l 16 (x, y) = xy(x -y) 2 , l 17 (x, y) = xy(x + y) 3 , l 18 (x, y) = xy(x -y) 3 , l 19 (x, y) = x 2 y 2 (x 2 + y 2 ). (6)
The 19 coefficients a i are obtained as follows. We consider the discrete values depending on ψ located at the eight points M k , 1 ≤ k ≤ 8, around the point M 0 , (see Fig. 2). From the discrete data at these points we determine 19 values to be interpolated by P(x, y) and its derivatives:

                                       Γ 1 (ψ) = ψ(M 1 ), Γ 2 (ψ) = ψ(M 2 ), Γ 3 (ψ) = ψ(M 3 ), Γ 4 (ψ) = ψ(M 4 ), Γ 5 (ψ) = ψ(M 0 ), Γ 6 (ψ) = ψ(M 5 ), Γ 7 (ψ) = ψ(M 6 ), Γ 8 (ψ) = ψ(M 7 ), Γ 9 (ψ) = ψ(M 8 ), Γ 10 (ψ) = (-∂ x + ∂ y ) ψ(M 1 ), Γ 11 (ψ) = ∂ y ψ(M 2 ), Γ 12 (ψ) = (∂ x + ∂ y ) ψ(M 3 ), Γ 13 (ψ) = -∂ x ψ(M 4 ), Γ 14 (ψ) = ∂ x ψ(M 0 ), Γ 15 (ψ) = ∂ y ψ(M 0 ), Γ 16 (ψ) = ∂ x ψ(M 5 ), Γ 17 (ψ) = (-∂ x -∂ y ) ψ(M 6 ), Γ 18 (ψ) = -∂ y ψ(M 7 ), Γ 19 (ψ) = (∂ x -∂ y ) ψ(M 8 ). (7)
Note that the derivatives at any point are taken in the direction of M 0 except that the full gradient is given at the point M 0 .

There is a one-to-one correspondence between the polynomial (5) and the above set of 19 data. More explicitely, the 19 coefficients a i in [START_REF] Chen | A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow[END_REF] are uniquely determined by the data (7). For the proof of this linear algebraic fact, see [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF].

In [START_REF] Chen | A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow[END_REF], the coefficients a i depend linearly on the data Γ k (ψ), 1 ≤ k ≤ 19. Therefore, P(x, y) can be rewritten as

P(x, y) = 19 ∑ i=1 19 ∑ j=1 A i, j Γ j (ψ) l i (x, y). (8)
We need to calculate the geometric coefficients A i, j , 1 ≤ i, j ≤ 19 in terms of the

vector h = [h 1 , h 2 , h 3 , h 4 , h 5 , h 6 , h 7 , h 8 ].
For this purpose, it is useful to decompose the polynomial P(x, y) into the sum of four terms P(x, y) = P(0, 0) + P 1 (x) + P 2 (y) + xyQ(x, y).

Looking at ( 5) and ( 6), these four terms are expressed as:

a 1 = P(0, 0) = ψ(M 0 ) ( given value) , (10)
   P 1 (x) = a 2 x + a 3 x 2 + a 4 x 3 + a 5 x 4 + a 6 x 5 ,
P 2 (y) = a 7 y + a 8 y 2 + a 9 y 3 + a 10 y 4 + a 11 y 5 .

(

) 11
The polynomial Q(x, y) in ( 9) is then defined as Q(x, y) = P(x, y) -P(0, 0) -P 1 (x) -P 2 (y) xy = a 12 + a 13 (x + y) + a 14 (xy) + a 15 (x + y) 2 + a 16 (xy) 2 (12)

+ a 17 (x + y) 3 + a 18 (x -y) 3 + a 19 xy(x 2 + y 2 ). ( 13
)
This decompsition is directional in the following sense:

• The polynomial P 1 (x) ∈ Span{x, x 2 , x 3 , x 4 , x 5 } corresponds to the "horizontal data". It is determined by the 5 data (see Fig. 2):

ψ(M 4 ), ψ(M 5 ), ∂ x ψ(M 4 ), ∂ x ψ(M 0 ), ∂ x ψ(M 5 ). (14)
• Similarly, P 2 (y) ∈ Span{y, y 2 , y 3 , y 4 , y 5 } corresponds to the "vertical data". It is specified by the 5 data

ψ(M 7 ), ψ(M 2 ), ∂ y ψ(M 7 ), ∂ y ψ(M 0 ), ∂ y ψ(M 2 ). (15)
• Finally, it can be shown that the polynomial Q(x, y) is determined by the 8 "diagonal data" in (7). These data are:

       Γ 1 (ψ) = ψ(M 1 ), Γ 3 (ψ) = ψ(M 3 ), Γ 7 (ψ) = ψ(M 6 ), Γ 9 (ψ) = ψ(M 8 ), Γ 10 (ψ) = (-∂ x + ∂ y )ψ(M 1 ), Γ 12 (ψ) = (∂ x + ∂ y)ψ(M 3 ), Γ 17 (ψ) = (-∂ x -∂ y )ψ(M 6 ), Γ 19 (ψ) = (∂ x -∂ y )ψ(M 8 ). ( 16
)
4 The numerical scheme

The embedded discrete biharmonic operator

In this section, we assume given for each point of the Cartesian grid the polynomial P(x, y) (5) in terms of the data Γ k (ψ). As explained in Section 3, the polynomial P M 0 (x, y) in ( 5) is explicitely known by the coefficients a i , given as the analytical functions:

h, [Γ j (ψ)] j=1,...,19 → a = a 1 , a 2 , . . . , a 18 , a 19 T . (17)
The discrete biharmonic at M 0 (x 0 , y 0 ) is obtained by:

∆ 2 h ψ(M 0 ) = 19 ∑ k=1 a k ∆ 2 l k (x 0 , y 0 ). (18)
There are four nonvanishing terms in the right-hand-side of (18) which are:

∆ 2 l 5 (x 0 , y 0 ) = 24, ∆ 2 l 10 (x 0 , y 0 ) = 24, ∆ 2 l 15 (x 0 , y 0 ) = 16, ∆ 2 l 16 (x 0 , y 0 ) = -16. ( 19
)
Therefore the discrete biharmonic at M 0 is given in terms of the coefficients

a k h, [Γ j (ψ)] by ∆ 2 h ψ(M 0 ) 24 a 5 (h, [Γ k (ψ)]) + a 10 (h, [Γ k (ψ)]) (20)
+ 16 a 15 (h, [Γ k (ψ)]) -a 16 (h, [Γ k (ψ)]) . (21)
The discrete equation at point M 0 is therefore (see (4)):

∆ 2 h ψ(M 0 ) = f (M 0 ). (22)
Equation ( 22) has to be supplemented by some additional relation connecting the derivatives ψ x,i, j , ψ y,i, j and the values ψ i, j . Our choice [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF][START_REF] Ben-Artzi | Navier-Stokes equations in planar domains[END_REF] is to use an Hermitian relation in the xand the ydirection. In the xdirection we have:

α 1,i ψ x,i-1, j + ψ x,i, j + α 2,i ψ x,i+1, j = β 1,i ψ i-1, j + β 2,i ψ i, j + β 3,i ψ i+1, j . (23)
The five coefficients α 1,i , α 2,i , β 1,i , β 2,i and β 3,i are defined as follows. Let M 0 = Q i, j (x i , y j ) and let the two neighboor points M 4 and M 5 be (see Fig. 2):

M 4 (x i -h i , y j ), M 5 (x i + h i+1 , y j ). ( 24
) Then      α 1,i = h 2 i+1 (h i+1 +h i ) 2 , α 2,i = h 2 i (h i+1 +h i ) 2 , β 2,i = 2h 4 i+1 +4h 3 i+1 h i -4h i+1 h 3 i -2h 4 i h i+1 (h i+1 +h i ) 3 h i , β 1,i = - 2h 4 i+1 +4h 3 i+1 h i h i+1 (h i+1 +h i ) 3 h i , β 3,i = 2h 4 i +4h i+1 h 3 i h i+1 (h i+1 +h i ) 3 h i . ( 25
)
In the ydirection we have γ 1, j ψ y,i, j-1 + ψ y,i, j + γ 2, j ψ y,i, j+1 = δ 1, j ψ i, j-1 + δ 2, j ψ i, j + δ 3, j ψ i, j+1 .

with values of the five coefficients γ 1, j , γ 2, j , δ 1, j , δ 2, j and δ 3, j deduced from the points M 7 and M 2 in a way similar to (25). We refer to [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF][START_REF] Ben-Artzi | Navier-Stokes equations in planar domains[END_REF] for an analysis of the Hermitian relations ( 23) and (26).

Assembling the global linear system

To each point (i, j) corresponds the discrete biharmonic relation ( 22) together with the horizontal and vertical Hermitian relations for the discrete gradient ( 23) and ( 26). All these relations form a linear system

AΨ = F. (27)
Assembling the matrix A using the relations (22, 23,26) is analogous to assembling the global matrix in the finite element method. According to Section 1, each point M i, j of the Cartesian grid belongs to one of the five categories:

1. interior regular calculated point 2. interior irregular calculated point 3. interior edge point 4. boundary point 5. exterior point In our computation, this classification is performed using a so-called level set model for the boundary ∂ Ω . Assume that (x, y) → ϕ(x, y) is a smooth function such that, at least locally ϕ(x, y)

     < 0 if (x, y) ∈ Ω , (interior point), > 0 if (x, y) ∈ Ω c , (exterior point), = 0 if (x, y) ∈ ∂ Ω , (boundary point).
(28) Following [START_REF] Chen | A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow[END_REF], the interior point

M 0 = M i, j is declared close to ∂ Ω if ϕ min,i, j ϕ max,i, j < 0 where ϕ min,i, j = min(ϕ i-1, j , ϕ i+1, j , ϕ i, j+1 , ϕ i, j-1 , ϕ i, j ), ϕ max,i, j = max(ϕ i-1, j , ϕ i+1, j , ϕ i, j+1 , ϕ i, j-1 , ϕ i, j ). (29)
In this case, the following quadratic model for ϕ is defined around M 0 by:

ϕ(x) = ϕ 0 + (∇ϕ 0 ) T .(x -x 0 ) + 1 2 (x -x 0 ) T (D 2 ϕ 0 )(x -x 0 ). (30)
In (30), ∇ϕ 0 and D 2 ϕ 0 stand for approximate values of the gradient and the Hessian of ϕ(x) at M 0 . In the computations, centered differences for ∇ϕ 0 and D 2 ϕ 0 are used. Using the model (30) allows to determine the approximate projection M * 0 of the interior point M 0 on ∂ Ω , [START_REF] Chen | A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow[END_REF]. This gives

M 0 = calculated point if dist(M 0 , M * 0 ) ≥ ε edge , edge point if dist(M 0 , M * 0 ) < ε edge . (31)
where ε edge is a fixed parameter. For each calculated point M 0 , the length vector h ∈ R 8 and the elementary matrix A i, j (h) ∈ M 19 (R) are evaluated according to the preceding classification into regular/irregular calculated points. Finally the elements of each matrix A i, j (h) are collected in the global matrix A. In a second step, for each point M i, j , the submatrix of A corresponding to the Hermitian relations for the derivatives ψ x and ψ y in ( 23) is calculated. The global linear system Aψ = b is the discrete version of the problem [START_REF] Abarbanel | Asymptotically Stable Fourth-Order Accurate Schemes for the Diffusion Equation on Complex Shapes[END_REF]. Note that it is solved by a direct solver. Fast solvers issues in the fashion of [START_REF] Chen | A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow[END_REF][START_REF] Ben-Artzi | A fast direct solver for the biharmonic problem in a rectangular grid[END_REF] will be addressed in a future work.

Numerical results

We present several numerical results for the biharmonic problem with additional Laplacian term:

α∆ 2 ψ(x) -β ∆ ψ(x) = f , x ∈ Ω , ψ = g 1 (x), ∂ ψ ∂ n = g 2 (x), x ∈ ∂ Ω . (32)
In each case, the domain Ω and the solution ψ(x) are specified. The right-hand side f (x) and the two boundary functions g 1 (x) and g 2 (x) are determined accordingly. The numerical scheme is then used to obtain an approximation for ψ based on the discrete values of f ,

Test cases in an ellipse

We first consider two test cases where the computational domain is an ellipse. A similar test case has already been considered in [START_REF] Ben-Artzi | A compact difference scheme for the biharmonic equation in planar irregular domains[END_REF]. The observed accuracy is very good. The order of convergence is located approximately in the interval I = [START_REF] Ben-Artzi | A fast direct solver for the biharmonic problem in a rectangular grid[END_REF][START_REF] Ben-Artzi | Navier-Stokes equations in planar domains[END_REF]. 1 Compact scheme for ∆ 2 ψ = f . The solution is ψ(x, y) = (1x 2 ) 2 (1y 2 ) 2 in the ellipse x 2 /1 2 + y 2 /2 2 ≤ 1. The ellipse parameters are (a = 1, b = 2, r = 1). The ellipse is embedded in the square [-2, 2] × [-2, 2]. We present e and e x , the l 2 errors for the streamfunction and for ∂ x ψ. The parameter for points labelled as edge points is ε edge = 5.10 We present e and e x , the l 2 errors for the streamfunction and for ∂ x ψ. The parameter for points labelled as edge points is ε edge = 5.10 -3 h.

Test cases in non convex domains

Star shaped domains

We consider first the biharmonic problem (see Example 4.3 in [5])

∆ 2 ψ(x) = 0 x ∈ Ω , ψ = g 1 (x), ∂ ψ ∂ n = g 2 (x), x ∈ ∂ Ω . (33)
The boundary of the domain is given in polar coordinates by

x(θ ) = R(θ ) cos(θ ), y(θ ) = R(θ ) sin(θ ), 0 ≤ θ < 2π, (34)
with R(θ ) = 0.6 + 0.25 sin(k p θ ). The domain is represented on Fig. 3 for k p = 7, (seven branches case). The exact solution is ψ(x, y) = x 2 + y 2 + e x cos(y). The numerical results are reported on Fig. 4 where the least square slope is represented, based on six grids. They show excellent accuracy, even for very coarse grids. Observe in addition the low error level for ψ and ∂ x ψ. The numerical results are reported on Fig. 6. Again, the accuracy is very good. But the levels of error are higher than in the flower case. This can be attributed to the non regular boundary.

Fig. 2

 2 Fig. 2 Zoom on the neighborhood of point M 0 on Fig. 1. The coordinates have been moved such that M 0 is the coordinates center. The 8 neighbors points of M 0 are the points M 1 , M 2 , M 3 , M 4 , M 5 , M 6 , M 7 and M 8 . The points M 1 , M 4 , M 6 and M 7 belong to the Cartesian grid. The points M 2 , M 3 , M 5 and M 8 belong to the boundary of the domain. They are obtained as the intersection of rays emanating from M 0 and directed towards M 2 , M 3 , M 5 and M 8 respectively. The points M 3 , M 5 and M 8 are outside the domain. The edge point above M 0 is marked with an open circle.

  )ψ = f . The solution is ψ(x, y) = 100(x 3 ln(1 + y)) + y 1+x in the ellipse (x -0.5) 2 /(0.5) 2 + (y -0.5) 2 /0.3 2 ≤ 1. The ellipse parameters are (a = 0.5, b = 0.3, r = 1) with center (x c , y c ) = (0.5, 0.5). The ellipse is embedded in the square [0, 1] × [0, 1].

Fig. 3

 3 Fig. 3 Seven branches star shaped domain embedded in a 33 × 33 grid. • Left: domain and grid. The boundary points are marked with black triangles. The edge points are marked with open circles. • Right: approximate solution corresponding to ψ ex (x, y) = x 2 + y 2 + e x cos(y).

Fig. 4

 4 Fig. 4 Star shaped domain: linear regression of the convergence rate for Ψ -(ψ ex ) ∞ and Ψ x -(ψ x,ex ) ∞ where the exact solution is ψ ex (x, y) = x 2 + y 2 + e x cos(y). • Left: domain with 7 branches, (k p = 7). • Right: domain with 9 branches, (k p = 9). On each regression line, the six points correspond to the six grids 10 × 10, 20 × 20, 30 × 30, 40 × 40, 50 × 50 and 60 × 60.

Fig. 5

 5 Fig. 5 Double circle shaped domain embedded in a 41 × 41 grid. Left: domain and grid. The boundary points are marked with black triangles. The edge points are marked with open cicles. Right: approximate solution ψ(x, y) = exp(x + y).

Fig. 6

 6 Fig. 6 Double circle shaped domain: linear regression and convergence rate for Ψ -(ψ ex ) ∞ and Ψ x -(ψ x,ex ) ∞ with • Left: ψ ex (x, y) = exp(x+y) • Right: ψ ex (x, y) = 10(x 5 sin(4πy)+ y 4 1+x 2 ). For each regression line, the six points correspond to the six grids 10 × 10, 20 × 20, 30 × 30, 40 × 40, 50 × 50 and 60 × 60.

  -3 h.

	mesh	17 × 17 Rate 33 × 33 Rate 65 × 65 Rate 129 × 129
	e ∞ 6.9555(-6) 3.53 6.0000(-7) 4.43 2.7790(-8) 3.29 2.8334(-9)
	(e x ) ∞ 4.0042(-4) 3.64 3.2033(-5) 4.07 1.9102(-6) 2.98 2.4215(-7)
	e 2 1.1759(-6) 3.15 1.3240(-7) 4.26 6.9034(-9) 3.44 6.3715(-10)
	(e x ) 2 7.4850(-5) 3.79 5.3933(-6) 3.98 3.4163(-7) 3.90 2.2865 (-8)

Table 2

 2 Compact scheme for ( 1 2 ∆ -∆ 2