Abe Junko

Abe Junko Reprentation

Junko Abe

REPRENTATION OF TAMASHII IN REFLECTIONS ON CONTEMPORARY TALES

Keywords: tamashii (spirit or soul), contemporary tales, vision of the death Matsutani

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

tales 日本民話の会. The stories she collected include several genres--tales, legends, myths, and so-called contemporary tales 現代の民話, which are the focus of this paper.

Matsutani gathered contemporary tales in quite an original way. In 1978, she founded the magazine Cahier of tales『民話の手帖』, in which a prepaid postcard was inserted. Readers were invited to write down a story they knew and send it to the Cahier of tales.

Each issue of the magazine had a theme. The first was the tengu 天狗, a supernatural creature with wings that lives in the mountains. The idea was a success, as people sent hundreds of stories, which would constitute the basis of the twelve volumes of Reflections on contemporary tales (1985-1996) 3 .

Each volume of the series contains from one to three different themes. We first have the traditional recurrent themes of Japanese tales, such as the yôkai 妖怪 (tengu, kappa 河童, wolf 狼, tanuki 狸, a sort of raccoon, cat 猫), and the kami kakushi (kidnapping by kami). Two volumes were dedicated to the Pacific war, because a lot of stories are related to this tragic moment of Japanese history. Three volumes deal with mysterious things or entities occurring in modern devices or inventions such as photography, television, trains, recording studios and so forth. The two final volumes are about tamashii 魂 tales.

Let us first give a definition of tamashii. In Japanese, tamashii corresponds to "spirit" or "soul", but this differs from the English term.

First, the tamashii resides in the body of a human being.

Second, unlike in the West where the soul does not leave a living body, the tamashii does. In some tales, the person is conscious of the departure of his or her tamashii and reports what he or she sees or hears while the tamashii is flying away. The tamashii that left the body can be seen in the form of a will-o'-the-wisp 火の玉 or the person him/herself.

The tamashii of a dead person takes the form of a will-o'-the-wisp or of a ghost. This interpretation is close to that of a ghost or phantom in Western culture.

This study is concerned with several chapters dedicated to the tamashii. I believe that the tamashii stories provide a unique opportunity to observe and better understand contemporary beliefs and ideas surrounding life and death in Japanese popular culture, which exist outside of the institutional religions, i.e. Buddhism and

Shintô.

The following chapters of Reflections on contemporary tales were chosen as the primary body for this research: "Will-o'-the-wisp"「火の玉」(volume 4), "Out-of-body experience「抜け出した魂」" (vol.4), "Announcement of death"「死の知らせ」 (vol.5), "Visit to the hereafter"「あの世へ行った話」(vol.5), "Supernatural stories from the frontline"「軍隊の怪談」(vol.2) .

The term "contemporary tales 現代の民話" could be misleading, since a "tale" is generally assumed to be something old, which has been transmitted orally from generation to generation. Matsutani Miyoko comments on the "contemporary tales" as follows:

These are stories told like this, as they come, but following well-defined patterns, as tales did in ancient times. They spread within one region, and people say "these things really happened."4 So these stories are definitely contemporary but firmly rooted in tradition.

Initially, the stories start as rumours, then turn into tales that follow traditional patterns even when their contents are modern. For example, at the beginning of the Meiji era, after trains were introduced to Japan in 1872, stories concerning false trains began to appear. The legendary tanuki or kitsune, who usually transform themselves, were said to have changed into trains running at full speed and to have collided against a real train and crashed. Shortly after, this story was transformed and became the history of transformation contest (bakekurabe 化けくらべ)5 .

Historically, Matsutani Miyoko is not the first writer to use the term "contemporary tale". The first mention is by the dramatist Kinoshita Junji 木 下 順 二 (1914-2006) who used the term in an essay called "Minwa kanken"「民話管見」(My humble thoughts on the folk tales), published in the journal Bungaku 『 文 学 』 (Literature) in 1952 6 . Kinoshita thought that he could use the idea for the framework of a play he was writing about contemporary society. Since everybody knew the story he wanted to use, he assumed that his play would acquire the strength of tradition.

Shortly after Kinoshita, another writer, Yamashiro Tomoe 山 代 巴 (1912-2004), started collecting "contemporary tales." 7 We can also trace back the idea of "contemporary tale" to 1910, when the Tales of Tôno 『 遠 野 物 語 』8 by Yanagita Kunio 柳 田 國 男 (1875-1962), were published.

Matsutani Miyoko integrated the Tales of Tôno in her Reflections on contemporary tales, claiming that they were nothing but contemporary tales.

The will-o'-the-wisp is not considered as the incarnation of the soul in the West, cf. for instance the novel by Drieu La Rochelle named The Will-o'-the-wisp (Le Feu follet). In The Devil's Pool (La Mare au diable) by George Sand, it is merely a natural phenomenon. On the contrary, in Japan, the will-o'-the-wisp is a physical manifestation of a tamashii. Even if it is commonly called hi no tama (fire ball) in standard language, various other names can be found in each region, attesting to its presence throughout Japan: tamashii, tamashi, hitodama, aonoro, hikarimono, chudama, chuntamashi, tamagai, shônendama, hikaridama, kechibi, kitsunebi or onibi. The morphology of will-o'-the-wisps is quite characteristic: they are red, orange, yellow, light blue or brown, and they fly away followed by a tail like shooting stars.9 I shall now analyse the representation of tamashii in Reflections of contemporary tales. I have identified six types of stories.

1) The appearance of tamashii, the harbinger of death

In different cultures, people have recounted strange phenomena heralding someone's imminent death. In The Legend of Death (La Légende de la mort) by Anatole Le Braz (1859-1926), a collection of Breton tales of the 19th century, a chapter10 is devoted to the so-called intersigns (in French intersignes, the Breton word being seblanchou or sinalieu): a bowl, a plate, a glass falls down by itself and breaks; a rooster crows in the afternoon; a magpie stands on the roof of a house; a dog howls to death during the night. These are some of the numerous signs for the forthcoming death of a relative, friend or neighbour.

Matsutani has also collected such stories in which death is announced by the same kind of intersigns: strangely croaking ravens; a Buddhist altar bell sounding by itself; a kakemono falling from the wall. 11 In both cultures these are signs of death. The person to whom such signs appear is rarely the one who will die.

However, the cases that interest us here are the appearance of the tamashii of the person who is going to die. This can be either through the will-o'-the-wisp or the shape of the person.

The most notable story involving a will-o'-the-wisp is "The Will-o'-the-wisp of Shipwrecked sinners" by IWASAKI Toshie of Miyagi Prefecture:

It must have occurred during the 35 th year of the Meiji era. A stranger rented the field in front of the Sasunohama cemetery. He mowed some stubble in and built a hut to live in, to cook and work. People report that every night, he used to hear the voices of many people, coming from the cemetery, going "ho, ho, ho, ho". They also spoke: "blah, blah, blah, blah." Every night… -What is that? wondered the man, going out of his cabin and looking towards the cemetery. Then he saw bright things flying around above the tombs: hyoi, hyoi, hyoi, hyoi.

-How strange, it really looks like a lot of people! It was a group of will-o'-the-wisp hovering around in a hubbub of voices and screams. The man then announced to the villagers of Sasunohama that a misfortune would happen. Shortly afterwards, the village's fishing vessel was wrecked on a reef in Arahama, and almost all men of the village disappeared.

After their funeral, the will-o'-the-wisp never showed up again. You see: it was indeed a premonition. 12 All stories like this follow a common scheme: someone sees the will-o'-the-wisp, they know that this is a harbinger of bad news and that something fatal is going to occur.

Invariably, at the end, one or more people die.

The storyteller tells this story in a lively way, using onomatopea and mimetic words.

A similar story is "Procession at the Ten.yû Temple " also related by IWASAKI.

In 1946, in Ogatsu village, a procession of will-o'-the-wisp used to go to the cemetery of the Ten.yû Temple every night. The villagers wondered what would happen. Shortly thereafter, a boat shuttling between surrounding villages was wrecked and thirty-eight villagers died. After the funeral, the will-o'-the-wisp never returned. 13 There are many kinds of stories involving will-o'-the-wisps in Reflections on contemporary tales: the will-o'-the-wisp leaves the house of a dying person and goes to the cemetery; or it stays outside a house or appears to relatives or friends from a great distance... The tamashii also may take the shape of a person. Testimonies abound of families of soldiers who died in the Pacific War, seeing their son, brother or grandson, appear at the moment of death, which was called gyokusai 玉 砕 (literally, a gem that breaks into a thousand pieces, an expression coined by the authorities to beautify their death).

2) The departing spirit says good-bye to its loved ones and friends

In the second type of story, the tamashii comes out of the body just before death to say goodbye to its loved ones and friends. They enjoy time together, and shortly afterwards, it becomes clear died at the very same time.

In this case, compared to the first type of stories, the tamashii has a good time with his or her family or friends. It is more than the simple announcement of death. This appears clearly in the following story, told by a resident of Kanagawa Prefecture:

Grandmother told me this story. It occurred around the 30th year of the Showa era. Not far from our house lived a sick old woman, who had been confined to bed for a long time. Every day, they helped her to take a bath, it was her great pleasure. One evening, we heard in our house: "Good evening, I came to take my bath." Clearly, this was the voice of the old woman, however, she was quite unable to walk up to our house… Grandmother listened carefully without moving from her bed. It was late. She heard the sound of the bath-tub lid taken away. A little later, it was put back. Then these words. "Thank you very much for the excellent bath! " Finally grandmother went up and saw the old woman leaving, dressed in a white yukata. She disappeared gradually. Two days later, people heard that the old woman had died. TAKAHASHI Hideko, Kanagawa prefecture 14 Note that in Japan, someone may invite you to their house to take a bath just as they would for a meal.

In one of the Tales of Tôno, a tofu maker named Masa tells his father's death story. He had been away for a long time because of a serious illness, and could not move at all. He suddenly appeared on the communal works of Shimotochinai where he used to work. He greeted everyone, worked in the fields with others and came back home at dusk. At night, the workers heard that he had died at the time he was with them. They realized that they had been working near his tamashii. 15 In numerous stories, the tamashii of the dying person goes to the family temple.

For example, in the Tales of Tôno, a man goes to see the monk of the temple. The monk is surprised to see him, but receives him with courtesy and offers him tea. They talk, then the man leaves. Finding the behaviour of his visitor a little odd, the monk asks one of his young disciples to follow him. The man is seen walking in the street, greeting several people, then turning at the corner of a street, where he disappears. His death was announced that evening. The monk later noticed that the tea had fallen between the tatami joint where the man had put his cup. 16 Stories of soldiers' tamashii spending time with their mother at the time of their death are particularly poignant.

Toshio, who fought at the front, appeared suddenly to his mother while she was sitting in the living room. He came to lay his head on her lap and fell asleep. His mother took the lice from his hair. Then he left. It was at that moment that he died at the Battle of 5 th May in 1945 at Okinawa. 17 Another soldier came home reciting the sutras in front of the house's altar.

Another one arrived hungry. When his mother turned back with the rice balls and a cup of green tea that she had prepared for him, he had gone.

Yet another soldier arrived in his military uniform, completely soaked. At that time he died on a sinking boat. 18 An eminent anthropologist, Iwata Shigenori (1961-), wrote an essay about this phenomenon of tamashii leaving the body just before death. He explained that the departure of tamashii from a living body only happens when the person is dying. The announcement of the death signifies the reunification of the tamashii and the body after spending time with relatives or friends. Iwata Shigenori also remarks that in Japan there exist rituals after the funeral that involve breaking rice bowls or putting upside down a mortar called usu 臼, so that the tamashii cannot lodge in the hollow of one of these containers and is thus able to reunite with the body. 19 When we hear these stories, we think that such phenomena cannot exist. To me, this is irrelevant. The point is that the protagonists believe that they really happened and that they hold such a concept of the tamashii.

3) The tamashii of a dead person comes to pick up a departing tamashii In the third type of story, the tamashii of a dead person comes to pick up a departing tamashii to lead it to the hereafter. A person in his or her deathbed starts talking to a deceased person. This may be their father, grandmother or a close friend.

Obviously, the two prepare his or her imminent departure. Soon afterwards, the person dies.

Reflections on contemporary tales contains the following story, related by KURIHARA Tadatoshi of Kagoshima prefecture.

Grandmother's health got worse, and all the family gathered in her house.

My father and I were in a room next to her's. It was raining and snowing at the same time. In the afternoon, suddenly, we heard the steps of a person coming barefoot to the house. However, we saw nobody outside. the dying often called the tamashii by his or her name, showing real intimacy: "Mom! ", "Grandfather!", "I am arriving!" There exists a real link between the two, and the atmosphere is always very warm.

Neither Buddha nor kami appear in this sort of tale. No representation of the hereafter does either, be it the Christian Paradise, the Yomi no kuni (the underground world of the deads in ancient Japan) or the Buddhist Pure Land of the West.

Besides the stories of omukae by a deceased family member, there are a few in which people (hito, onna no hito, rojin) come down from the sky, although these are rather ambiguous. We do not know exactly who these people are, nor where they come from.

Sometimes the family members who are attending the scene, understanding what is going to happen, proceed to the ritual of the call-of-tamashii (tamayobai 魂呼ば い). They shout loudly the soon-to-decease person's name. Some add: "Don't go now!", "It's too early!".

In some stories, a horse carriage, a boat or even a car picks up the dying person.

In other cases, will-o'-the-wisp can embody the tamashii of the deceased. There is no conversation. The will-o'-the-wisp comes down from a sacred mountain or the cemetery of the family temple. It enters the house through the chimney and shortly after, two will-o'-the-wisp return to the sacred mountain. People say, for example: the grandfather came to pick up the grandmother.

So we do not die alone. Our relatives protect us, they keep close to us, they assist us. This makes death less terrifying.

4) Near-death experience

The forth type of story is the near-death experience. Somebody close to death finds him or herself in the middle of a field of flowers. There, they come across several dead people who beg him not to approach them. At that very moment, the person hears people shouting his or her name and comes back to life.

Tales of Tôno, collected in Reflections on the contemporary tales, is one example. Near-death experiences usually occur when people who have been clinically dead for a few moments return to life because of medicine. In the United States, this phenomenon has been observed countless times by specialists, who have collected many testimonials. Almost all follow the same pattern, though the details differ according to the cultures.

At the village of

In Japan, especially in stories collected by Matsutani Miyoko, the images are usually of fields of poppies, sunflowers or other flowers, in which are many dead people.

Between our world and the hereafter, a river almost always marks the frontier. The dead are on the other side. Sometimes they say do not come. In other cases, they try to take you to the other side of the river.

"I am glad to learn that my father and my brother are in such a nice place. If so, it doesn't bother me to go there", a woman said after such an experience.

5) The reincarnation

In Reflections on contemporary tales, two types of reincarnations can be identified. First, we encounter the type in which people write calligraphy of a sutra or the name of the deceased on the palm of his or her hand or on the sole of his or her foot.

Sometimes, a baby is born in the village with the same inscription on its hand or foot.

The family says that the baby is the reincarnation of the dead person. The inscription disappears when the new born takes its first bath in which water some soil taken from the grave of the deceased has been added.

About ten stories of this type are collected in Reflections on contemporary tales.

Here is an example.

At the end of the Meiji era, in the village of Yamagata, there was a beggar named Okume, who often paid visit to an important family of the village named K. When Okume died, the family called the bonze to organize her funeral. She thus got a grave. For the peace of her tamashii, the bonze wrote a sutra on the sole of her foot.

Later, a girl was born in the family. She had the same inscription on her foot. They bathed the baby adding in the water a small amount of soil taken from Okume's tomb. The sutra disappeared at once. Thereafter, the daughter had a happy destiny.

Shiobara Tsuneko, Nagano prefecture 24

The anthropologist Miyata Noboru (1936-2000) says that in Japan, when someone dies, it is customary after the funeral to wish for his or her reincarnation.

During the Mikokiki ミコキキ, a death ritual which aims at bringing down the tamashii of the deceased into the body of the miko 巫女 (medium), the relatives of the deceased person are happy when they hear his or her voice say, "I will be reborn soon in the person of a little child." Moreover, in popular traditions, when a grandfather or a grandmother dies, and that a birth occurs immediately afterwards in the family, people consider that the deceased person has been reincarnated. 25 When an elderly person dies, it seems that it is not customary to write sutra, but in any case, if a child is born in the family, people rejoice themselves because they think that it is the reincarnation of the grandmother or grandfather.

The second type of reincarnation stories has to do with children who are supposed to remember their previous life, being able to provide details about it, mentioning people or places names.

So did a woman from the Tokyo prefecture, who remembered her childhood in the 1970s. "I suddenly started to say that I had been Masako. I could remember very clearly my life as Masako. "Masako was the sister of her father, she had died at the age of five. Her father was bewildered as she told things that only Masako could have told. 26 In another story, a doctor realized, during a hypnosis session of a girl of thirteen years old named Junko, that in her subconscience, she had precise remembrance of the life of another little girl from Tokyo, Kikuko, who had died very young. The tamashii of Kikuko had left her body just before the incineration, had flied away to enter inside a calf named Jirobei at Hiroshima. After, this tamashii came to Okayama in the family of Orihara to be reborn as Junko. 27 This type of story is known through the Chronicle of the Rebirth of Katsugorô 『勝五郎の再生記聞』 (1828) by Hirata Atsutane 平田篤胤 (1776-1843), a philosopher of National Studies 国学 of the Edo period. Katsugorô, who lived at Hachiôji in Tokyo, told at the age of 8 to his parents that he had been the child of a peasant family. This boy had died prematurely. Katsugorô remembered that his tamashii had remained first on a desk at home, and then, led by an old man with white hair, had been reincarnated.

Lafcadio Hearn (1850-1904), wrote also about Katsugorô in the short story "Rebirth of Katsugorô"「勝五郎の再生」in 1897.

This concept of reincarnation is different from the Buddhist idea of transmigration. In Buddhism, behaviour during one's life has influence on the reincarnation. Thus, in case of misbehaviour, one reincarnates as an insect or an animal, sign of punishment. As we saw earlier according to the anthropologist Miyata Noboru's writings, reincarnation follows, in Japanese popular traditions, the will of the family who would like the deceased to be reborn. 28 This idea of reincarnation is reassuring.

The deceased continues to exist through the life of another person.

6) The Wondering or errant tamashii

Finally, the sixth and final type of stories comprises those about wandering or errant souls. In this case, the tamashii appears only after death.

Persons who died tragically, violently from an accident, illness, suicide or in the war, appears to living people, because their tamashii do not find peace 浮かばれない.

their homes, fearing they might be drawn into the procession. 30 In Japan, even today, if people see this kind of phenomenon, they proceed to a ceremony in order to appease the tamashii. They invite a monk who will recite the sutra, or they erect a stupa.

One Japanese specificity resides in stories where suddenly deceased persons feel regrets about unaccomplished acts during their life.

In April 1993, the wife of the writer Takashi Yoichi died suddenly. 1999. 31 Op.cit, vol.5, pp.417-419 In such stories, the dead person's regrets disappear when a living person understands that the spirit of the dead is still around, and comforts the spirit: asking a monk to chant a sutra, organizing a memorial service or building a tomb or a memorial tower. So, the end of this story is an exception.

★ ★ ★

In Japanese, we can say "to go on a journey", to mean "die." "My grandfather died," therefore, gives "My grandfather went on a journey 祖父は旅立ちました"。 Symbolically, it means that death is not only an end but also a starting point.

The popular Japanese stories show how the living and the dead are linked in close relationships, emotionally and even physically.

After death, as the body disappears, nothing is left. This is the materialistic vision of death according to modern rationalism. It generally results in fear, and sometimes in an extreme feeling of solitude. In comparison, the image of death in Japanese popular traditions is very different.

The Portuguese writer Wenceslau de Moraes (1854-1929), who lived at the beginning of the 20 th -century on the island of Shikoku, became very fond of this concept of the well-meaning, guiding omnipresence of tamashii. He admired the Japanese, for whom the deceased persons stay near their family and loved ones. 32 In Europe, a similar concept of death disappeared long time ago.

32 DE MORAES Wenceslau, O « Bon-odori » em Tokushima『徳島の盆踊り』, translated from portuguese by OKAMURA Takiko, Kôdansha, Gakujutsu bunko, 1998, pp.179-186

  In this story, only the grandmother can hear the voice of Kesa don. The rest of the family only hears the footsteps of the uncle who stopped in front of the grandmother's

	grandmother's room.
	"Ah Kesa-don! Have you come to get me? Thank you."
	Grandmother said.
	"You know, I am not ready yet. You can sleep here tonight. And we will go
	away tomorrow morning. "
	Kesa-don is the name of grandmother's dead brother.
	Grandmother goes on chatting in a low voice.
	She died early in the morning. We said to each other: "Our uncle came to find
	grandmother. " 20
	the other could not see." 21 This demonstrates the persistence of popular beliefs
	concerning the proximity of the afterlife and the close back-and-forth connections
	between the two worlds.
	In about twenty stories of omukae from Reflections on the contemporary tales 22 ,
	The footsteps reached

room, but not his voice. We also note the relaxed tone of the grandmother, who sounds happy to go on a trip with the brother she loved so much. She greets him as if she had left him the day before. "It's been a long time! Have you come? "... This phenomenon is commonly called omukae お迎え, the picking up of a dying person to guide him or her to the hereafter. Apart from contemporary tales, we also have some recent testimony. A survey conducted by a team of young doctors, nurses, philosophers, and psychologists working with elderly people and targeting 350 families revealed that 42 percent of the subjects interviewed testified that just before his or her death, one of their relatives "talked about the existence of persons or landscapes that

20 Op.cit, vol.5,

SHIMIZU Tetsurô (editor) 清水哲郎監修、How to live, how to die, Thanatology from the hospice 『どう生き どう死ぬか、現場から考える死生学』, Kyûsen

Shoin, 2009, p.163. 22 Op.cit, vol.5, pp.130-149.

  One month later, Yoichi saw one of these large black butterflies called swallowtail butterfly flying around him. At another moment, he noticed that the same butterfly was sitting on the unfinished painting of his wife representing the Kirishima mountain. For Yoichi, this was the sign of a feeling of regret that his wife was communicating to him. He decided to bring back the painting to the bottom of the Kirishima mountain, where they had a house. Just as he sat in his car, he saw the black butterfly sitting on his wife's sit ! He started. After a few kilometers, the butterfly fluttered in all directions, as if it wanted him to go faster. His wife loved speed, and she used to ask him to drive faster. Upon arrival, Yoichi placed the painting in front of the mountain. The butterfly followed him and went on the painting. Then it disappeared. When Yoichi took again the wheel, the butterfly was dead on the seat of the passenger. He made for it a tomb and implored it not to return never again. It happened five or six years ago, but the butterfly returned more than once.Takashi Yoichi, interviewed by Matsutani Miyoko, Fukuoka prefecture 31

30

Itô Susumu 伊藤進、The forest and the Devil『森と悪魔』 、Iwanami shoten, 2002, p.245.

Cf. also Chasse fantastiques et cohortes de la nuit au moyen âge, Claude Lecouteux, Imago,

Reflections on contemporary tales『現代民話考』, 12 volumes, Chikuma bunko, Chikuma shobô, 1985-1996.

Contemporary tales『現代の民話』, MATSUTANI Miyoko, Chûkô shinsho, Chûo kôron, p.15.

Ibid.,

"Minwa Kanken「民話管見」(My humble thoughts on the folk tales)", Bungaku『文学』 (Literature), may 1952, Iwanami shoten.

YAMASHIRO Tomoe 山代巴、Creators of tales『民話を生む人々』, Iwanami Shoten, 1958.

YANAGITA Kunio 柳田國男、Tales ofTôno『遠野物語』, 1910.

Reflections on contemporary tales, vol.4,

The Legend of Death (La Légende de la mort), Anatole Le Braz, tome1, 1912, Champion, pp.1-68.

Op.cit. vol.5,

MATSUTANI Miyoko 松谷みよ子・編著、Tales of Onagawa and Ogatsu, recounted by IWASAKI Toshie 『女川・雄勝の民話・岩崎としゑの語り』,Kokudosha, 1987, pp.99-102.

Ibid.,

Tale No. 86, Op.cit.

Tale No. 87, Op.cit.

Op.cit, vol.5,

Ibid,

IWATA Shigenori 岩田重則, "Reflections on the tamashii that leaves the body"「遊離霊魂論」 in KOMATSU Kazuhiko 小松和彦 (editor), Views of the parallel world in Japan『日本人の異界 観』 、Serika shobô,

2006, p.188.

Tale No. 97, Op.cit, 1910.

MIYATA Noboru 宮田登、Folklore of spirit and trip『霊魂と旅のフォークロア』, 宮田登 日 本を語る Vol. 7、YoshikawaKôbunkan, 2006, pp.78-80.

Op.cit, vol. 5, p525.

Ibid., p.531.

Op.cit, vol.2, p.402.

They are full of regret and bitterness.

Many stories are related during the Pacific War, when young people were massively killed. Here is a example of stories of dead soldiers who came back as ghosts.

On the island of Guadalcanal, during a ground offensive by 15,000 US marines, the Ikki regiment lost its 800 fighters only in nine days. Fifteen days later, the regiment Kawaguchi arrived in the island. At night, two Japanese sentries suddenly heard the noise of boots on the sand. They saw fifteen soldiers coming toward them from the seaside. Holding the Japanese flag, their head down, injured, limping heavily, the troop passed before them and disappeared into the darkness. When the sentries told this story, their chief did not believe them. So he joined them the next day and they saw the same scene. It was the ghosts of the Ikki regiment.

Narrated by Masuda Katsumi, Tokyo 29

These stories of errant soldiers' tamashii occur many times in Reflections on contemporary tales.

We can make an interesting parallel between these Japanese stories and legends of western ghosts and phantoms. In literature, let us remember the famous scene in the beginning of Shakespeare's play Hamlet, where the phantom of Hamlet's father appears. Horatio, Hamlet's friend, verifies the apparition that the two guards, Bernardo and Marcellus, said they had seen on previous nights. The phantom passes once more near these three men (Act I, Scene I). We know that the father had been tragically killed by his brother. His soul could not rest in peace and wandered endlessly.

We can also cite the "ghost army" of the Middle Ages. Dead soldiers, led by a skeleton, cross the sky at great speed. The peasants used to confine themselves inside


Carine David

Victor David

Victor David Management

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

MANAGEMENT REGIMES

Carine David, University of New Caledonia

Victor David,IRD Nouméa New Caledonia is an autonomous local authority within the Republic of France, with its own territorial government, congress and provincial assemblies. The political "Noumea Accord" of 1998 initiated a staged transfer of most of administrative powers from France to New Caledonia.

Importantly, the 3 provinces (South Province, North Province and Loyalty Islands Province) have jurisdiction over all matters that are not expressly vested in the State or territory. As environment is not mentioned within the French State or New Caledonia's power, most of environmental regulations are adopted by the provincial assemblies.

Nevertheless, some important powers regarding the law of the sea are vested in New Caledonia, directly or indirectly related to protection and conservation of marine species, areas or habitats, including public lands, customary land ownership; public health, zoosanitary and phytosanitary safety, exploitation and management of natural resources in the exclusive economic zone (including deep sea mining), traffic and road transport and ports.

The French State has primary jurisdiction for international affairs and as such most of international conventions ratified by France are applicable in New Caledonia (i.e. Montego Bay Convention on the law of the sea of 1982). The autonomous territory can also sign and ratify regional conventions. Importantly, in July 2008, the World Heritage Commission listed 'The lagoons of New Caledonia: reef diversity and associated ecosystems' as a serial listing of 6 marine clusters containing 15743 km² of coral reefs.

In 2014, the government of New Caledonia created the 'Natural Park of the Coral Sea' with the ambition of protecting the ecological integrity of the marine environment. The main objectives of the Natural Park of the Coral Sea are the protection of ecosystems and the maintaining of the existing links between them, preservation of rare and migratory endangered species and identification and enhancement of the Natural Parks cultural heritage. In 2017, a management plan of the Natural Park of the Coral Sea was presented to the government, in order to achieve these objectives.

The division of powers between the different levels of governance has resulted in a great complexity regarding the responsibilities of each level of government in relation to environmental matters. Moreover, another important feature is the strong customary tenure and management practices of the Kanak people on marine areas, even if the ownership of these are vested within the provinces.

Regarding protection and conservation measures and legal protective instruments related to marine species, areas and habitats, most of the regulations and the management are however the responsibility of provincial authorities.

Each province adopted its own Code of environment, respectively in 2008 (North province [START_REF]Délibération n°306-2008[END_REF] ), 2009 (South province [START_REF]3 Délibération 2016-13/API 2016 portant adoption du[END_REF] ) and 2016 (Loyalty Islands Province 3 ).

The structure of the three codes is very similar, but there are significant differences in their substantive provisions. In each Code, there are regulations dealing with the preservation of natural heritage, including protected areas, ecosystems and species. There are also provisions about natural resource management, including fisheries, and, in the case of South Province, exploitation of biological, genetic and biochemical resources 4 .

Protected Areas

The protected area provisions of North and South provinces code were substantially updated when the codes were adopted. In both provinces, protected areas have been consolidated into several categories, aligned with the internationally recognised IUCN protected area categories [START_REF]Guidelines for Applying Protected Area Management Categories[END_REF] . Each type of protected areas has clearly defined objectives and management rules. These protected area types and existing protected areas have been reclassified into four categories in the South Province [START_REF]Strict nature reserve (IUCN Category Ia), provincial park (IUCN Category II), nature reserve (IUCN Category IV) and sustainable resource management area (IUCN Category VI)[END_REF] and seven categories in the North province [START_REF]Strict nature reserve (IUCN Category Ia), wilderness reserve (IUCN Category 1b), provincial park (IUCN Category II), natural or cultural heritage site (IUCN Category III), nature reserves (IUCN Category IV), natural or cultural landscape (IUCN Category V), sustainable resource management area (IUCN Category VI)[END_REF] . The protected areas regulation which will soon be enacted in the Loyalty Islands province, on the contrary, is not based on IUCN categories but is seeking to formalise customary practices within marine areas traditionally managed by kanak population.

Ecosystem Protection

The South Province code provides protection for some 'ecosystems of heritage interest'. Among these, mangroves, seagrass beds and coral reefs are listed ecosystems. Regarding such ecosystems, the regulation requires provincial authorisation in any case an operation is likely to have a significant impact on them. Such an authorisation is only issued if an impact study demonstrates that the operation will not significantly undermine the conservation status of the ecosystem. However, in the absence of an alternative, authorisation may be granted on 'public interest' grounds, subject to mitigation or compensation conditions. As the notion of 'public interest' is not defined, many authorisations are granted on this ground.

Neither the North Province, nor the Loyalty Islands province codes provide for such protection of ecosystems.

Species Protection

There are, in North and South provinces, existing regulations protecting endemic, rare and threatened species are in force, in order to respond the severe threats to New Caledonia's biodiversity. Listed protected species include whales, dolphins, dugongs, marine turtles, birds, fishes, freshwater crustaceans, molluscs.

It is prohibited to hunt, capture, destroy, remove, disturb, transport, use, trade or possess a protected animal species, or to destroy or remove its eggs or nest. It is also prohibited to destroy, modify, alter or degrade the habitat of a protected species. Provincial authorisation may be issued for collecting, possessing or trading a protected species for conservation or educational purposes, including breeding, research or presentation to the public.

The Loyalty Islands province has not yet enacted such a regulation, but some regulations protect a few species, such as coconut crabs.

The provincial codes do not deal with international trade in endangered species, as implementation of the Convention on International Trade in Endangered Species (CITES) as there is a legal debate whether this matter falls within the jurisdiction of the territory or the provinces. As a result, the Convention has not yet been implemented, even if the French State administration is taking great care about the respect of the stipulations of the Convention as it is its responsibility that can be engaged.

Invasive Species

Invasive species can have an important impact on the New Caledonian environment. Each code provides for a list of invasive plant and animal species. The codes prohibit wilful, negligent or reckless introduction of listed invasive species into the wild. It is also prohibited to cultivate, transport, use, sell or purchase any invasive plant. The provincial administration may authorise the production, use or possession of listed invasive species for agricultural, fisheries or forestry purposes (crayfish for example) or, in the Loyalty Islands province, for customary practices.

Fishing 8

The fisheries regulations in North and South provinces have been substantially modified during the codification process. In each province, the code emphasises the 'longterm preservation' and 'responsible and rational exploitation' of fisheries resources. Professional and artisanal fishers are required to apply for a fishing licence, and recreational fishing is regulated. Fishing methods and gear are controlled, size limits apply and conservation measures are defined for certain marine species (e.g. groupers, lobsters, mud crabs, clams, trocas, oysters, bêche-de-mer). Further management measures may be established (e.g. total allowable catch, individual quotas, seasonal closures, area closures, gear restrictions). North Province specifically recognises customary marine resource management practices and undertakes to work with customary resource owners to more fully integrate these practices in future regulations. The Loyalty Islands province is planning to enact a new regulation with the same concerns about customary practices.

Clark P, David C., New provincial environmental legislation in New Caledonia : continuity and reform in environmental governance in a French pacific territory », Asia Pacific Journal of Environmental Law, n° 13, 2011, pp. 135-145


Jonas Koko Limos

Efficient MATLAB codes for the 2D/3D Stokes equation with the mini-element

Keywords: Finite element method, Stokes problem, Uzawa conjugate gradient, MAT-LAB AMS subject classification: 65N30, 76D07, 65K10

We propose a fast MATLAB implementation of the mini-element (i.e. P 1-Bubble/P 1) for the finite element approximation of the generalized Stokes equation in 2D and 3D. We use cell arrays to derive vectorized assembling functions. We also propose a Uzawa conjugate gradient method as an iterative solver for the global Stokes system. Numerical experiments show that our implementation has an (almost) optimal time-scaling. For 3D problems, the proposed Uzawa conjugate gradient algorithm outperforms MAT-LAB built-in linear solvers.

Introduction

MATLAB is a popular problem solving environment, widely used for general scientific computation in education, engineering and research. MATLAB is nowadays a standard tools in many areas. Thanks to its collection of direct (e.g. LU , LDL , Cholesky) and iterative (e.g. conjugate gradient, GMRES, bi-conjugate gradient) solvers, it is tempting to use MATLAB for the numerical approximation of partial differential equations (PDEs). For the finite element method (FEM), the main obstacle for using MATLAB is the assembly of the global matrix and vector. Since MATLAB built-in solvers are optimized, the assembly operations can take up to 99% of whole CPU time, as shown in [START_REF] Koko | Vectorized MATLAB codes for two-dimensional linear elasticity[END_REF], using an implementation with standard loop over triangles ( [START_REF] Alberty | Remarks around 50 lines of matlab: short finite element implementation[END_REF][START_REF] Alberty | Matlab implementation of the finite element method in elasticity[END_REF][START_REF] Kwon | The Finite Element Method Using MATLAB[END_REF]) directly derived from compiled languages (FORTRAN and C/C++). Unfortunately, some vectorized FEM codes are less flexible and require a huge amount of memory due to the allocation of auxiliary arrays and the corresponding matrix operations [START_REF] Koko | Vectorized MATLAB codes for two-dimensional linear elasticity[END_REF][START_REF] Rahman | Fast MATLAB assembly of FEM matrices in 2D and 3D: Nodal elements[END_REF][START_REF] Rahman | Fast MATLAB assembly of FEM matrices in 2D and 3D: Edge elements[END_REF]. Recently Koko [START_REF] Koko | Fast MATLAB assembly of fem matrices in 2d and 3d using cell array approach[END_REF] proposes a MATLAB 1 implementation close to the standard form by using cell-arrays to store the gradient of the basis functions, for the Poisson equation and linear elasticity in 2D and 3D.

In this paper, we propose a fast MATLAB implementation of the P 1-Bubble/P 1 finite element (Mini element, [START_REF] Arnold | A stable finite element for the Stokes equations[END_REF][START_REF] Ern | Eléments finis: théorie, applications, mise en oeuvre[END_REF][START_REF] Glowinski | Numerical methods for fluids (part 3)[END_REF]) for the generalized Stokes problem in 2D and 3D. The mini element for spatial discretization of the Stokes problem is easy to use in engineering practice since its allows for the use of equal-order interpolation (the same mesh for velocity and pressure). Equal-order interpolation is very useful in large-scale multi-physics codes. Indeed, a code dealing with several independent variables (e.g. chemical species, velocity components, etc) requires the transfer of information between its different components at each time-step. Fast implementation means that our code operates on array and does not use for-loops over elements (triangles of tetrahedrons) for the assembling operations. Instead, we use cell-arrays to store element matrices as in [START_REF] Koko | Fast MATLAB assembly of fem matrices in 2d and 3d using cell array approach[END_REF]. We also propose a solution strategy for the final linear system. Indeed, we propose a efficient (preconditioned) Uzawa conjugate gradient method derived from the one used with P 2/P 1 (or P 1-iso-P 2/P 1) finite element pair ( [START_REF] Cahouet | Some fast 3-D solvers for the generalized Stokes problem[END_REF][START_REF] Glowinski | Numerical methods for fluids (part 3)[END_REF]). The proposed conjugate gradient method operates on the dual (pressure) space and, at each iteration, d independent linear systems are solved (d = 2, 3). Our implementation needs only MATLAB basic distribution functions and can be easily modified and refined.

The paper is organized as follows. The model problem is described in Section 2, followed by a finite element discretization in Section 3. The element matrices in 2D/3D are described in Section 4. In Section 5, we propose our Uzawa conjugate gradient method for solving the Stokes system. MATLAB implementation details are given in Section 6. Numerical experiments are carried out in Section 7. Readers can download and edit the codes from http://www.isima.fr/~jkoko/Codes/KSTOK.tar.gz.

The model problem

Let Ω be a bounded domain in R d (d = 2, 3) with a Lipschitz-continuous boundary Γ. Consider in Ω the Stokes problem 3) turn to be a generalized Stokes problem encountered in time discretization of Navier-Stokes equations (see e.g. [START_REF] Glowinski | Numerical methods for fluids (part 3)[END_REF]). The constant ν > 0 is the kinematic viscosity.

αu -ν∆u + ∇p = f , in Ω, ( 2.1)
∇ • u = 0, in Ω, (2.2)
u = u D , on Γ, (2.3
We need the functional spaces

V = H 1 0 (Ω) d , V D = v ∈ H 1 (Ω) d : v = u D on Γ , P = q ∈ L 2 (Ω) : Ω q dx = 0 ,
and bilinear forms

a i (u i , v i ) = α(u i , v i ) Ω + ν(∇u i , ∇v i ) Ω , i = 1, . . . , d a(u, v) = d i=1 a i (u i , v i ),
where (•, •) Ω stands for the standard L 2 (Ω) scalar product. The variational formulation of the Stokes problem (2.1)-(2.3) is as follows:

Find (u, p) ∈ V D × P such that:

a(u, v) -(p, ∇ • v) Ω = (f , v) Ω , ∀v ∈ V , (2.4)
-(q, ∇ • u) Ω = 0, ∀q ∈ P.

(2.5)

3 Finite element discretization with P 1-Bubble-P 1

For the finite element discretization of (2.4)-(2.5), we have to chose a finite element pair for the velocity field and the pressure. This choice cannot be arbitrary but must satisfy the inf-sup condition [START_REF] Babuska | Error bounds for finite element method[END_REF][START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers[END_REF].

Mini element

In this paper we study the discretization of the Stokes problem (2.1)-(2.3) by the finite element pair P 1-bubble/P 1 (the so-called mini-element), introduced by Arnold, Brezzi and Fortin [START_REF] Arnold | A stable finite element for the Stokes equations[END_REF]. This element leads to a relatively low number of degrees of freedom with a good approximate solution.

Let T h be a triangulation of Ω and T a triangle of T h . We define the space associated with the bubble by

B h = v h ∈ C 0 ( Ω); ∀T ∈ T h , v h|T = xb (T ) .
We also defined the discrete function spaces

V ih = v h ∈ C 0 ( Ω); v h|T ∈ P 1 , ∀T ∈ T h ; v h|Γ = 0 , i = 1, . . . , d P h = q h ∈ C 0 ( Ω); q h|T ∈ P 1 , ∀T ∈ T h ; : Ω q h dx = 0 ,
and we set

X ih = V ih ⊕ B h and X h = X 1h × • • • × X dh .
With the above preparations, the discrete variational problem reads as follows.

Find (u h , p h ) ∈ X h × P h such that a(u h , v h ) -(p h , ∇ • v h ) Ω = (f , v h ) Ω ∀v h ∈ X h , (3.1)
-(q h , ∇ • u h ) Ω = 0 ∀q h ∈ P h . (3.2)
For a given triangle T , the velocity field u h and the pressure p h are approximated by linear combinations of the basis functions in the form

u h (x) = d+1 i=1 φ i (x)u i + u b φ b (x), p h (x) = d+1 i=1 φ i (x)p i ,
where u i and p i are nodal values of u h and p h while u b is the bubble value. The basis functions are defined by

φ 1 (x) = 1 -x -y, φ 2 (x) = x, φ 3 (x) = y, φ b (x) = 27 3 i=1 φ i (x),
in 2D, and

φ 1 (x) = 1 -x -y -z, φ 2 (x) = x, φ 3 (x) = y, φ 4 (x) = z, φ b (x) = 256 4 i=1 φ i (x) in 3D.
To construct the coefficient matrices, a number of integrals involving powers of the basis functions will be computed. Integrals over a triangle T can be evaluated directly by the following formulas

T φ α 1 1 φ α 2 2 φ α 3 3 dx = 2|T | α 1 ! α 2 ! α 3 ! (α 1 + α 2 + α 3 + 2)! (3.3) T φ α 1 1 φ α 2 2 φ α 3 3 φ α 4 4 dx = 6|T | α 1 ! α 2 ! α 3 ! α 4 (α 1 + α 2 + α 3 + α 4 + 3)! , (3.4)
where |T | stands for the triangle area (in 2D), or the tetrahedron volume (in 3D). A useful property for the basis functions is

d+1 i=1 ∇φ i = 0. (3.5)

Algebraic formulation

We use the following notations for the discrete velocity/pressure nodal values ūi =

u i u ib , fi = f i f ib , i = 1, . . . , d. (3.6)
System (3.1)-(3.2) leads to the following algebraic form, using notations (3.6),

   Ā 0 -B 1 0 Ā -B 2 -B1 -B2 0       ū1 ū2 p    =    f1 f2 0    (3.7) in 2D, or       Ā 0 0 -B 1 0 Ā 0 -B 2 0 0 Ā -B 3 -B1 -B2 -B3 0            ū1 ū2 ū3 p      =      f1 f2 f3 0      (3.8)
in 3D. In (3.7) and (3.8) Ā = M + R, with M the mass matrix and R the stiffness matrix. Bi is the divergence submatrix associated with the i-th partial derivative, i.e.

Bi ≡ (q h , ∂ i u ih ) Ω , i = 1, . . . , d.

To create the algebraic system (3.7) or (3.8), the discrete system (3.1)-(3.2) is evaluated over each triangle T to obtain the element matrices and vectors

M (T ) ij = T αφ i φ j dx, R(T ) ij = T ν∇φ i • ∇φ j dx, B(T ) ij = T ∂ 1 φ i φ j dx + T ∂ 2 φ i φ j dx, f (T ) i = T f φ i dx.
Then assembling operations consist of direct-summing the element matrices over the triangulation T h to obtain the global matrices M = ( Mij ), R = ( Rij ), B = ( Bij ) and f = ( fi )

Mij = T ∈T h M (T ) ij , Rij = T ∈T h R(T ) ij , Bij = T ∈T h B(T ) ij , fi = T ∈T h f (T ) i .
In the next sections we detail the element matrices M , R, Bi and the right-hand side f .

Element matrices

For P 1 finite element, we need only the element area and the gradient of the basis functions for the element matrices and vectors. To derive vectorized MATLAB codes, we need analytical expressions forall element matrices and vectors. This section is devoted to this task. For the sake of the presentation we drop the (T )-superscript introduced in the previous section to distinguish element matrices from global ones.

Two-dimensional case

For a triangle T , let {(x i , y i )} i=1,2,3 be the vertices and {φ} i=1,2,3 the corresponding basis functions. The gradient of φ i are given by

   ∇φ t 1 ∇φ t 2 ∇φ t 3    =    1 1 1 x 1 x 2 x 3 y 1 y 2 y 3    -1   0 0 1 0 0 1   = 1 2|T |   y 2 -y 3 x 3 -x 2 y 3 -y 1 x 1 -x 3 y 1 -y 2 x 2 -x 1   , (4.1)
where |T | is the area of T given by

2|T | = det x 2 -x 1 x 3 -x 1 y 2 -y 1 y 3 -y 1 = (x 2 -x 1 )(y 3 -y 1 ) -(x 3 -x 1 )(y 2 -y 1 ).
A nonbubble entry of the element stiffness matrix R, is given by

Rij = T ν∇φ i • ∇φ j dx = |T |ν∇φ i • ∇φ j , i, j = 1, 2, 3. (4.2)
For the bubble entries Rbj , for j = 1, 2, 3. A straightforward calculation yields to (using (3.5))

Rbj = 9 4 |T | 3 i=1 ∇φ i = 0, j = 1, 2, 3.
For the diagonal entry corresponding to the bubble (i.e. i = j = b) we have

Rbb = ν T 27 2 ∇(φ 1 φ 2 φ 3 ) • ∇(φ 1 φ 2 φ 3 )dx = 81 10 ν|T | |∇φ 1 | 2 + |∇φ 2 | 2 + |∇φ 3 | 2 + ∇φ 1 ∇φ 2 + ∇φ 1 ∇φ 3 + ∇φ 2 ∇φ 3 = 81 10 ν|T | |∇φ 1 | 2 + |∇φ 2 | 2 + ∇φ 1 • ∇φ 2 =: ω R , (4.3)
using (3.5).

With the above results, the element stiffness matrix is therefore

R = R 0 0 ω R .
As for the stiffness matrix, we set M = ( Mij ) i,j=1,...,3 the nonbubble part of the mass matrix. A direct calculation yields

M ij = α 6 |T | if i = j, α 12 |T | elsewhere. (4.4)
The bubble part of the mass matrix is given by

Mbj = 3α 20 |T |, j = 1, 2, 3, Mbb = 81α 280 |T | =: ω M . (4.5)
The element mass matrix is therefore

M = M z z ω M , where z = 3 20 α|T | (1 1 1) .
Finally, the element stiffness/mass matrix Ā is

Ā = A z z ω .
where we have set

A = R + M and ω = ω R + ω M . A direct integration yields the element divergence matrix -B = [-B1 -B2 ], where Bi = [B i -B ib ] with B i = |T | 3    ∂ i φ 1 ∂ i φ 2 ∂ i φ 3 ∂ i φ 1 ∂ i φ 2 ∂ i φ 3 ∂ i φ 1 ∂ i φ 2 ∂ i φ 3    , i = 1, 2, (4.6)
and

B ib = 9|T | 20    ∂ i φ 1 ∂ i φ 2 ∂ i φ 3    , i = 1, 2. (4.7)
The contribution of the right-hand side component f i , in nonbubble terms, is given by

f (T ) i = |T | 3 f iT , i = 1, 2, 3 (4.8)
where f iT is a mean value of f i on T . The bubble terms of the right-hand side are

f (T ) ib = T f i φ b dx = 9 20 |T |f iT , i = 1, 2. (4.9)
With the above element matrices and vectors the 11 × 11 element system corresponding to

(3.7) is         A z 0 0 -B 1 z ω 0 0 B 1b 0 0 A z -B 2 0 0 z ω B 2b -B 1 B 1b -B 2 B 2b 0                 u 1 u 1b u 2 u 2b p         =         f 1 f 1b f 2 f 2b 0         .
To reveal diagonal blocks, the system above can be rearranged as follows

         A 0 z 0 -B 1 0 A 0 z -B 2 z 0 ω 0 B 1b 0 z 0 ω B 2b -B 1 -B 2 B 1b B 2b 0                  u 1 u 2 u 1b u 2b p         =         f 1 f 2 f 1b f 2b 0         . (4.10)
We can now eliminate the bubble unknowns u 1b and u 2b since they correspond to diagonal blocks (the ω blocks). From (4.10) 3 and (4.10) 4 , we deduce that

u ib = (f ib -B ib p -z u i )/ω, i = 1, 2. ( 4
   A -ω -1 zz 0 -B 1 -ω -1 zB 1b 0 A -ω -1 zz -B 2 -ω -1 zB 2b -B 1 -ω -1 B 1b z -B 2 -ω -1 B 2b z -ω -1 (B 1b B 1b + B 2b B 2b )    (4.12)
and the right-hand side

   f 1 -ω -1 zf 1b f 2 -ω -1 zf 2b -ω -1 (B 1b f 1b + B 2b f 2b )    . (4.13)

Three-dimensional case

Let {x i = (x i , y i , z i )} i=1,...,4 be the vertices of a tetrahedron T and {φ i } i=1,...,4 the corresponding basis functions. The gradient ∇ x φ i over T are given by

      ∇φ 1 ∇φ 2 ∇φ 3 ∇φ 4       =      1 1 1 1 x 1 x 2 x 3 x 4 y 1 y 2 y 3 y 4 z 1 z 2 z 3 z 4      -1      0 0 0 1 0 0 0 1 0 0 0 1     
An alternative formula for computing ∇φ i is

∇ x φ i = J -1 ∇ ξ φ i (ξ).
where J is the Jacobean matrix of the mapping ξ → x(ξ), that is,

J =   x 2 -x 1 y 2 -y 1 z 2 -z 1 x 3 -x 1 y 3 -y 1 z 3 -z 1 x 4 -x 1 y 4 -y 1 z 4 -z 1   . (4.14)
The volume of tetrahedron T is given, from (4.14), by |T | = det(J)/6. As for 2D case, the nonbubble entries of the element stiffness matrix ar given by

Rbj = 9 4 |T | 3 i=1 ∇φ i = 0, j = 1, 2, 3, 4. (4.15)
while Rbj = 0, for all j = 1, . . . , 4. For the diagonal entry, using (3.4)-(3.5), we obtain

Rbb = ν T 256 2 ∇(φ 1 φ 2 φ 3 φ 4 ) • ∇(φ 1 φ 2 φ 3 φ 4 )dx = 8192 845 ν|T | 3 =1 |∇φ | 2 + ∇φ 1 • ∇φ 2 + ∇φ 1 • ∇φ 3 + ∇φ 2 • ∇φ 3 =: ω R .(4.16)
For the 3D mass matrix, a straightforward calculation with a linear tetrahedron shows that (for nonbubble entries)

Mij =    α T |T | 10 if i = j, α T |T | 20 if i = j
, where α T is a mean value of α on T . The bubble part of the mass matrix is given by

Mbj = 8α 105 |T |, j = 1, 2, 3, Mbb = 8192α 51975 |T | =: ω M . (4.17) 10
The element mass matrix is therefore

M = M z z ω M , where z = 8α 105 |T | 1 1 1 1 . If we set -B = [-B1 -B2 -B3 ] and Bi = [B i -B ib ]
, a straightforward calculation yields to

B i = |T | 4       ∂ i φ 1 ∂ i φ 2 ∂ i φ 3 ∂ i φ 4 ∂ i φ 1 ∂ i φ 2 ∂ i φ 3 ∂ i φ 4 ∂ i φ 1 ∂ i φ 2 ∂ i φ 3 ∂ i φ 4 ∂ i φ 1 ∂ i φ 2 ∂ i φ 3 ∂ i φ 4       , i = 1, . . . , 4; (4.18)
and

B ib = 32 105 |T |     ∂ i φ 1 , ∂ i φ 2 ∂ i φ 3 ∂ i φ 4     . (4.19)
The contribution of the right-hand side component f i , in nonbubble terms, is given by

f (T ) i = |T | 4 f iT , i = 1, 2 , 3, 4 (4.20)
where f iT is a mean value of f i on T . The bubble terms of the right-hand side are

f (T ) ib = T f i φ b dx = 32 105 |T |f iT , i = 1, 2, 3. (4.21)
As in 2D, we rearrange the Stokes system and we get

u ib = (f ib -B ib p -z u i )/ω, i = 1, . . . , 4. (4.22)
Substituting (4.22) into the rest of the system, we obtain a linear system whose matrix is

         A -ω -1 zz 0 0 -B 1 -ω -1 zB 1b 0 A -ω -1 zz 0 -B 2 -ω -1 zB 2b 0 0 A -ω -1 zz -B 3 -ω -1 zB 3b -B 1 -ω -1 B 1b z -B 2 -ω -1 B 2b z -B 3 -ω -1 B 3b z -ω -1 3 i=1 B ib B ib          (4.23)
and the right-hand side

       f 1 -ω -1 zf 1b f 2 -ω -1 zf 2b -ω -1 3 i=1 B ib B ib        . (4.24)

Uzawa conjugate gradient Algorithm

We propose in this section, a preconditioned conjugate gradient method for solving the Stokes system after the assembly of the element systems (4.12)-(4.13) and (4.23)-(4.24). For a finite element pair of the form P k+1 /P k (e.g. P 2/P 1 or P 1-iso-P 2/P 1), the preconditioner advocated by [START_REF] Cahouet | Some fast 3-D solvers for the generalized Stokes problem[END_REF] (see also [START_REF] Fortin | Augmented Lagrangian Methods: Application to the Numerical Solution of Boundary-Value Problems[END_REF][START_REF] Glowinski | Numerical methods for fluids (part 3)[END_REF][START_REF] Glowinski | Augmented Lagrangian and Operator-splitting Methods in Nonlinear Mechanics[END_REF]) is efficient and widely used. In the author knwoledge, there is not an equivalent preconditioner for the pair P 1-Bubble/P 1 (or P 1/P 1 with stabilization). The algorithm presented in his section is therefore an original contribution.

The Stokes system can be rewritten in a compact form

A -B -B -C u p = f f p (5.1)
where A is (symmetric) positive definite block diagonal matrix, C is positive semi-definite matrix.

Let us introduce the generalized Lagrangian function (5.3)

L (u, p) = 1 2 u Au -f u -p Bu - 1 2 p Cp -f p p. ( 5

Uzawa Conjuage gradient algorithm

To derive a dual (Uzawa) algorithm for (5.1), we assume that u = u(p) is the solution of Poisson equation

Au = f + B p (5.4)
that is, in the decomposed useful form

A i u i = f i + B i p, i = 1, . . . , d.
Then multiplying (5.4) by u and substituting the result in (5.2), we obtain

L (u(p), p) = - 1 2 u Au - 1 2 p Cp -f p p. (5.5)
If we introduce the dual functional

J * (p) = -min u L (u(p), p)
with u(p) solution of (5.4), the saddle-point problem (5.3) becomes Find p such that J * (p) ≤ J * (q), ∀q.

(5.6) From (5.5), J * is quadratic and coercive. From (5.4), we deduce that the mapping p → u(p) is linear and u(p + td) = u(p) + tw where w is the solution of the sensitivity problem

Aw = B d, (5.7)
or

A i w i = B i d, i = 1, . . . , d.
It follows that, the derivative of

J * is r := ∇J * (p) = Bu + Cp + f p . (5.8)
With a search direction d, we compute an optimal stepsize ρ by solving

∇J * (p + ρd) d = 0, that is ρ = -d (Bw + Cd)/(r d),
where w is the solution of the sensitivity equation (5.7). Since J * is a quadratic function, the best search direction is a conjugate gradient direction. As a consequence, the best algorithm for (5.6) is a conjugate gradient algorithm. At each iteration k, the Fletcher-Reeves conjugate gradient direction is given by

d k = r k+1 + β k d k β k = r k+1 2 r k -2 .
A dual (Uzawa) conjugate gradient algorithm for solving the saddle-point problem (5.1) is outlined in Algorithm 1. Theoretically, Algorithm 1 converges in at most n B = rank(B) iterations. Obviously, for large scale problems, n B is very large and it is preferable that convergence should be obtained in a number of iterations considerably less than n B .

Algorithm 1 Uzawa conjugate gradient algorithm for (5.1)

k = 0. p 0 is given 0.1 Initial solution A i u 0 i = f i -B i p 0 , i = 1, . . . , d r 0 = Bu 0 + Cp 0 + f p 0.2 Initial gradient and direction r 0 = Bu 0 + Cp 0 + f p d 0 = r 0 . k ≥ 0. While g k g k > ε(g 0 g 0 ) k.

Sensitivity and stepsize

A i w k i = B i d k , i = 1, . . . , d r k = Bw k + Cd k ρ k = (r k d k ) (d k r k ) k.2 Update p k+1 = p k -ρ k d k , u k+1 = u k -ρ k w k k.3 New gradient r k+1 = r k -ρ k r k , k.4 Conjugate gradient direction β k = (r k+1 r k+1 )(r k r k ) -1 d k+1 = r k+1 + β k d k .

Preconditioned Uzawa conjugate gradient algorithm

A practical implementation of a conjugate gradient method for solving (5.1) requires a preconditioner, that is, a suitable scalar product for computing ∇J * (p) instead of the standard one used in (5.8). The convergence properties of the conjugate gradient method for the generalized Stokes problem are deteriorated for large values of the ratio α/ν, see e.g. [START_REF] Glowinski | Numerical methods for fluids (part 3)[END_REF][START_REF] Glowinski | Augmented Lagrangian and Operator-splitting Methods in Nonlinear Mechanics[END_REF].

To derive a preconditioner for (5.1) following the idea of [START_REF] Cahouet | Some fast 3-D solvers for the generalized Stokes problem[END_REF][START_REF] Glowinski | On the predoncitioned conjugate gradient solution of a Stokes problem with Robin-type boundary conditions[END_REF], we need to simplify the continuous problem. We first notice that the equivalence between P 1-Bubble/P 1 element and the stabilized formulation has been proved [START_REF] Pierre | Regularization procedures of mixed finite element approximation of the Stokes problem[END_REF][START_REF] Pierre | Regularization procedures of mixed finite element approximation of the Stokes problem[END_REF][START_REF] Matsumoto | A relationship between stabilized FEM and bubble functions element stabilization method with orthogonal basis for incompressible flows[END_REF][START_REF] Pierre | Optimal selection of the bubble function in the stabilization of P1-P1 element for the Stokes problem[END_REF][START_REF] Baiocchi | Virtual bubbles and Galerkin-Least-Squares methods[END_REF]. Then, if we neglect the bubble contribution in the stiffness and divergence matrices, (5.1) can be expressed in the strong form as αu -ν∆u + ∇p = f (5.9)

∇ • u -∇ • (ν h ∇u) = 0, (5.10)
where ν h is an element dependent constant.

As in [START_REF] Glowinski | On the predoncitioned conjugate gradient solution of a Stokes problem with Robin-type boundary conditions[END_REF] we define the linear operator from L 2 (Ω) into L 2 (Ω) (neglecting the constant term in (5.8))

φ := Aq = ∇ • u q -∇ • (ν h ∇q) (5.11)
where u q is the solution of αu q -ν∆u q = -∇q.

(5.12)

Note that (5.12) is the (strong) sensitivity system. The idea behind preconditioning is to find a linear operator B such that BAq = q. Applying the divergence operator in (5.12), we obtain α∇ • u q -ν∇ • ∆u q = -∆q

or α∇ • u q -ν∆(∇ • u q ) = -∆q (5.13)
Using (5.11) in (5.13), we obtain

-∆q -α∇ • (ν h ∇q) + ν∆ (∇ • (ν h ∇q)) = αφ -ν∆φ. (5.14)
Then, in practice, at each step of the preconditioned conjugate gradient algorithm, the gradient of J * is computed as an approximate solution of the linear system

(K + αC + νKM -1 C)g = (αM + νK)r (5.15)
where K and M are the stiffness and the mass matrices, respectively. Let us introduce the mesh Reynolds number

Re h = α ν h 2
where h is the mesh size. Taking into account the CPU time and the storage requirement for computing the last term of the matrix involved in (5.15), we consider the following preconditioning system instead (K + αC)g = Hr.

(5.16)

where

H = νK + αdiag(M ) if Re h ≤ 1 νK + αI if Re h > 1 (5.17)
In (5.16) and (5.17), K and M are P 1 stiffness and mass matrix, respectively, and C the bubble matrix from Function 6.1-6.3.

From the theory of preconditioned conjugate gradient methods for the Stokes problem (see e.g. [START_REF] Glowinski | Numerical methods for fluids (part 3)[END_REF]) if Dirichlet conditions are imposed for the velocity, then for g in (5.14) we must impose ∂g/∂n = 0 (homogeneous Neumann boundary conditions). On the other hand, where a stress condition is prescribed for the fluid, we must impose g = 0 in (5.14).

Algorithm 2 Preconditioned Uzawa conjugate gradient algorithm for (5.1) k = 0. p 0 is given: Set P = K + αC, H given by (5.17).

Initial solution and residual

A i u 0 i = f i -B i p 0 , i = 1, . . . , d r 0 = Bu 0 + Cp 0 + f p 0.2 Initial gradient and direction Pg 0 = Hr 0 d 0 = g 0 . k ≥ 0. While g k r k > ε(g 0 r 0 ) k.

Sensitivity and stepsize

A i w k i = B i d k , i = 1, . . . , d r k = Bw k + Cd k ρ k = (r k d k ) (d k r k ) k.2 Update p k+1 = p k -ρ k d k , r k+1 = r k -ρ k r k , u k+1 = u k -ρ k w k k.3 New gradient P g k = H r k , g k+1 = g k -ρ k g k k.4 Conjugate gradient direction β k = (g k+1 r k+1 )(g k r k ) -1 d k+1 = g k+1 + β k d k .
Remark 5.1. The preconditioning system (5.15), that is, (αM + νK) -1 (K + αC + νKC)g = r, induces (over the discrete pressure space) the norm |g| P = g r.

Remark 5.2. If C = 0 in 5.15, we recover the preconditioning system of Cahouet and Chabat [START_REF] Cahouet | Some fast 3-D solvers for the generalized Stokes problem[END_REF] widely used for P 2/P 1 or P 1-iso-P 2/P 1 finite elements.

With the preparation given in the previous section, the preconditioned Uzawa conjugate gradient algorithm for solving the Stokes system (5.1) is detailed in Algorithm 2. Algorithm 2 inherits the main properties from the conjugate gradient algorithm of Cahouet-Chabat [START_REF] Cahouet | Some fast 3-D solvers for the generalized Stokes problem[END_REF]:

• At each step, we solve d independent linear systems for the velocity field. The matrices involved are sparse and can be factorized or preconditioned once and for all.

• The same algorithm is able to deal with 2D or 3D systems without any complication.

MATLAB implementation

We know detail the assembly of the Stokes systems (4.12)-(4.13) and (4.23)-(4.24). For the computational efficiency, the MATLAB codes must be vectorized (i.e. without long forloops). We then use arrays, cell-arrays, and MATLAB vectorized operators and functions

.*, ./, .^, sum, sparse

Mesh representation and KPDE package

We assume that the triangulation of Ω consists of np nodes and nt elements (triangles or tetrahedrons). We adopt the mesh representation by arrays used in [START_REF] Koko | Vectorized MATLAB codes for two-dimensional linear elasticity[END_REF][START_REF] Koko | Fast MATLAB assembly of fem matrices in 2d and 3d using cell array approach[END_REF][START_REF] Persson | A simple mesh generator in Matlab[END_REF]. The nodes coordinates are stored in an array p(1:np,1:2) (in 2D) or p(1:np,1:3) (in 3D). The element nodes are stores in an array t(1:nt,1:3) (in 2D) or t(1:nt,1:4) (in 3D). Dirichlet boundary conditions are provided by a list of nodes and the corresponding prescribed values.

As shown in [START_REF] Koko | Fast MATLAB assembly of fem matrices in 2d and 3d using cell array approach[END_REF] using cell-array in FEM allows to have implementation close to the standard form used in classical languages (C/C++, FORTRAN) FEM codes, while being efficient. The idea is to compute and store, for all triangles, the element matrix entry A (T ) ij in the cell-array At{i,j}. Then we use MATLAB function sparse to assemble A with finite small for loops for i=1:nd for j=1:nd A=A+sparse(t(:,i),t(:,j),At{i,j},np,np); end end where nd=3 (in 2D) or nd=4 (in 3D). This approach is used in KPDE package [START_REF] Koko | Fast MATLAB assembly of fem matrices in 2d and 3d using cell array approach[END_REF] for assembling matrices and vectors from Poisson and linear elasticity equations in 2D and 3D.

We need two key functions from KPDE package, kpde2dgphi and kpde3dgphi, which compute the gradient of the basis functions and the elements volume as follows

[ar,g]=kpde2dgphi(p,t); %2D [vol,g]=kpde3dgphi(p,t); %3D
g is 3-by-1 or 4-by-1 cell-array such that g{i}(:,k) is ∂ k φ i (x) for all elements.

2D case

Using the cell array g and the array ar, the nonbubble entry Rij , (4.2), of the element stiffness matrix is then given (for all triangles) by Rij=nu*ar.*sum(g{i}.*g{j},2);

The bubble diagonal entry ω R , ( The entries (i, j) of the element divergence matrix (4.6) are vetorized as follows MATLAB Functions 6.1-6.2 assemble the Stokes system. To make the assembling functions self-contained, all calculations are integrated except elements area and the gradient of basis functions which can be computed outside and passed as argument. Note that in 2D, the time for computing the triangles area and the gradient of the basis functions is not significant. Functions 6.1-6.2 can be called without the last two arguments without a significant computational overcost.

Instead of matrix (4.12) and vector (4.13), Functions 6.1-6.2 can return submatrices and subvectors of (4.12)-(4.13) if called with more than one output arguments. The statement [A,B,C]=kstok2dp1bmat(p,t,nu,alpha) w(s)=R'\(R\b(s)); But for large scale 3D problems, computing R'\(R\b(s)) requires a large amount of memory and can fail.

Numerical experiments

We now propose some numerical experiments to demonstrate the performances of our implementation. The computations have been carried out on a Dell Precision T3610 work station equipped with Intel Xeon 3.0GHz processor with 32GB RAM. The MATLAB version is 9 (R2016a).

Scalability

We first study the scalability of our MATLAB codes: We consider the discretization of a unit cube (0, 1) d (d = 2, 3) with a uniform mesh of size h, with nt triangles (or tetrahedrons) and np nodes. This initial mesh is successively uniformly refined to produce meshes of size h/2, h/4, h/8, etc. After each refinement the number of triangles is multiplied by 4 (2D) and the number of tetrahedraons by 8 (3D). Since the assembly process is essentially based on the number of elements, we expect that the time to assemble the matrices increases by approximately the same factor, i.e. 5 -6 in 2D and 8 -10 in 3D as observed in [START_REF] Koko | Fast MATLAB assembly of fem matrices in 2d and 3d using cell array approach[END_REF] for Poisson equation and linear elasticity. Table 1-2 show the assembly CPU times (in Seconds) for the Stokes system in 2D and 3D, respectively. We can notice an almost linear optimal time-scaling for our implementation.

Factorization versus Uzawa conjugate gradient

We now consider a Stokes flow in a driven cavity Ω = (0, 1) d with f = 0 in (2.1) and 2D: Γ 1 = (0, 1) × {1}, u = (1, 0) on Γ 1 , and u = 0 on ∂Ω \ Γ 1 ;

3D: Γ 1 = (0, 1) × (0, 1) × {1}, u = (1, 0, 0) on Γ 1 , and u = 0 on ∂Ω \ Γ 1 .

Ω is discretized by uniform meshes of size 1/16, 1/32, 1/64, 1/128, 1/256 and 1/512 in 2D, and 1/4, 1/8, 1/16, 1/32 and 1/64 in 3D. Since thre are powerful linear (direct) solvers in MATLAB, we first compare our conjugate gradient algorithm with

• \ (backslash) the standard MATLAB solver for general matrix based on Gaussian elimination;

• ldl the block LDL factorization for symmetric indefinite systems.

Table 3 shows the comparative performances for the 2D Stokes problem with ν = 1/50 and α/ν = 10 3 . We can notice that the proposed Uzawa conjugate gradient (Algorithm 2) and the MATLAB uilt-in Gaussian elimination are almost equivalent even though, in Algorithm 2, the component systems are uncoupled and can be solved in parallel. The CPU times for LDL include CPU timde for the factorization and columns and rows permutation to reduce fill-in that represents up to 90% of the whole CPU time.

If the Stokes problem is used in an iterative process (e.g. time stepping or linearization), then Algorithm 2 or LDL factorization are preferable. Indeed, if a LDL is carried out (once and for all) in the initialization step, then the solution of linear system reduces to forward/backward substitutions in the rest of the iterative process. The computational cost of Algorithm 2 can be reduced by using, as initial solution at the current step, the solution of the previous step. For the 3D Stokes problem the proposed Algorithm 2 outperforms the MATLAB direct solvers, Table 4. For the largest problem (4 * 65 3 = 1 098 500 unknowns) the direct solvers fail because of lack of memory due to fill-in during factorization. Table 5 shows the performances of the MATLAB gmres iterative solver using incomplete LU factorization as preconditioner (MATLAB function ilu with 10 -3 as drop tolerance). The value of restart parameter is 10. GMRES algorithm outperforms Algorithm 2 up to h = 1/32. But for the largest problem, Algorithm 2 is more than four times faster. It is clear that for large scale 3D problems, the proposed Uzawa algorithm is preferable.

2D Visualization

In two-dimensional incompressible fluid problems, it is usual to display the stream-lines. If the domain Ω is bounded and simply connected, in order to compute the stream-function ψ, we have to solve the Poisson-Neumann problem

-∆ψ = ω, in Ω, (7.1)
∂ n ψ = -u • τ, (7.2)
where ω = ∂ 1 u 2 -∂ 2 u 1 is the vorticity and τ the counter-clockwise oriented unit tangent vector at Γ. Problem (7.1)-( 7.2) has a unique solution in H 1 (Ω)/R. The variational formulation of (7.1)-(7.2) is Find psi ∈ H 1 (Ω) (∇ψ, ∇ϕ) Ω = (u 1 , ∂ 2 ϕ) Ω -(u 2 , ∂ 1 ϕ) Ω , ∀ϕ ∈ H 1 (Ω), (7.3) This leads to the following algebraic system using P 1 finite element

Rψ = B 2 u 1 -B 1 u 2
where R is the 2D Laplacian matrix. We impose ψ = 0 at an arbitrary node to ensure the uniqueness. We know consider a test problem derived from a benchmark problem described in [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF]. A mesh sample is shown in Figure 1. The inflow and outflow conditions (on left/right boundaries) are u 1 = 0.3 0.41 2 × 4y(0.41 -y), u 2 = 0 on Γ in = {0} × (0, 0.41), u 1 = 0.3 0.41 2 × 4y(0.41 -y), u 2 = 0 on Γ out = {2.2} × (0, 0.41).

On the other parts of the boundary of Ω, homogeneous boundary conditions are precribed (i.e. u = 0). The parameter α in (2.1) is set to 0. The center of the internal cylinder is (0.25, 0.2) and the diameter is 0.1. The kinematic viscosity is ν = 10 -3 . This gives a Reynolds number of Re = 30 based on the diameter of the cylinder and the maximum of the inflow velocity. The domain is discretized by a non uniform mesh consisting of 1730 nodes and 3280 triangles, Figure 1. The velocity field obtained with MATLAB command Unfortunately, for 3D flows, there is no simple tool for graphics output. quiver3 allows for visualisation of 3D velocity fields but the result is often unsatisfactory. Plotting 3D functions (or their contours) is a non trivial problem. There is no simple subproblem like (7.3) for streamlines in 3D.

Conclusion

We have proposed a fast MATLAB package for the numerical approximation of the generalized Stokes problem with the mini-element. Numerical experiments show that the proposed assembling functions have an optimal linear time-scaling. The proposed Uzawa conjugate gradient algorithm outperforms the MATLAB built-in solvers for 3D problems.

  ) where u = (u 1 , . . . , u d ) ∈ R d is the velocity vector, p the pressure and f = (f 1 , . . . , f d ) ∈ R d the field of external forces. In equation (2.1), α ≥ 0 is an arbitrary constant. If α = 0, then equations (2.1)-(2.3) turn to be the classic Stokes problem. If α > 0, then equations (2.1)-(2.

. 2 )

 2 Due to the properties of A and C, the saddle-point for (5.2) exists. Then it follows that (5.1) is the saddle-point equation for (5.2), that is, (5.1) characterizes the solution of the saddle-point problem min u max p L (u, p) = max p min u L (u, p).

  B1ij=ar.*g{i}(:,1)/3; B2ij=ar.*g{i}(:,2)/3; while the bubble entries (4.7) are B1ib=(9/10)*ar.*g{i}(:,1); B2ib=(9/10)*ar.*g{i}(:,2); The vectorized element contribution of the right-hand side (4.8) is f1=(1/3)*ar.*f1t; f2=(1/3)*ar.*f2t; while the bubble entries (4.9) are f1b=(9/20)*ar.*f1t; f2b=(9/20)*ar.*f2t;

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Mesh sample of the flow around cylinder problem

Figure 4 :

 4 Figure 4: Streamlines for the flow around a cylinder problem

Table 1 :

 1 CPU times (in Seconds) for assembling the Stokes matrix system of size N in 2D.

	h	1/32	1/64	1/128	1/256	1/512	1/1024
	N 3*1089 3*4225 3*16641 3*66049 3*263169 3*1050625
	A	0.018	0.071	0.325	1.516	7.242	36.390
	b	0.002	0.011	0.050	0.251	1.214	5.995
	h	1/4	1/8	1/16	1/32	1/64	1/128
	N 4*125 4*729 4*4913 4*35937 4*274625 4*2146689
	A 0.075 0.057	0.397	3.438	34.538	348.249
	b	0.027 0.018	0.051	0.465	4.544	46.394

Table 2 :

 2 CPU times (in Seconds) for assembling the Stokes matrix system of size N in 3D.

Table 3 :

 3 Comparative performances of MATLAB direct solvers and Algorithm 2 for the 2D Stokes system, ν = 1/50 and α/ν = 10 3

Table 4 :

 4 Table 6 shows the good convergence properties of the proposed Uzawa conjugate gradient algorithm when α/ν >> 1. Comparative CPU times (in Sec.) of MATLAB direct solvers and Algorithm 2 for the 3D Stokes system, ν = 1/50 and α/ν = 10 3 , (OoM= Out of Memory).

	Mesh size h	1/4 1/8 1/16 1/32	1/64
	Gaussian elim. 0.01 0.05 1.73 145.44 OoM
	LDL	0.02 0.08 6.29 128.54 OoM
	Algorithm 2	0.04 0.89 1.70 78.08	1894.59
	Mesh size h	1/4 1/8 1/16 1/32 1/64
	CPU Times (Sec.) 0.04 0.02 1.66 39.61 8230.59

Table 5 :

 5 Performances of the gmres iterative solver for the 3D Stokes system, ν = 1/50 and α/ν = 10 3 α/ν 10 1 10 2 10 3 10 4 10 5 10 6 10 7

	ν = 1/50	44	39	41	28	13	8	7
	ν = 1/200	37	34	33	26	13	8	7
	ν = 1/1000 31	28	27	24	18	8	7

Table 6 :

 6 Number of iterations versus α/ν for the 3D driven cavity problem for h = 1/32

Function 6.1 Assembly of the 2D Stokes matrix (4.12).

function [A,B,C]=kstok2dp1bmat(p,t,nu,alpha,ar,g) % A=kstok2dp1bmat(p,t,nu,alpha) or [A,B,C]=kstok2dp1bmat(p,t,nu,alpha) % A=kstok2dp1bmat(p,t,nu,alpha,ar,g) or [A,B,C]=kstok2dp1bmat(p,t,nu,alpha,ar,g) %- ---------------------------------------------------------------------------- returns the velocity right-hand side b=[b1; b2] and the pressure right-hand side bp. These submatrices are used in the preconditioned Uzawa conjugate gradient method presented in Section 5. A=[Ah Z Z Bh{1}'; Z Ah Z Bh{2}'; Z Z Ah Bh{3}'; Bh{1} Bh{2} Bh{3} Ch]; elseif (nargout == 3) A=Ah; B=[-Bh{1} -Bh{2} -Bh{3}]; C=-Ch; end MATLAB Functions 6.3-6.4 assemble the three-dimensional Stokes system. for computational efficiency, the elements volume vol and the gradient of the basis functions g can be computed once and for all, and passed to Functions 6.3-6.4. For more flexibility, vol and g can also be computed inside Functions 6.3-6.4 if they do not appear in the list of input arguments. As in 2D, Functions 6.3-6.4 can return submatrices and subvectors used in the preconditioned Uzawa conjugate gradient algorithm, if called with more than one output argument. Function 6.4 Assembly of the 3D Stokes right-hand side (4.24) function [b,bp]=kstok3dp1brhs(p,t,f1,f2,f3,nu,alpha,vol,g) % b=kstok3dp1brhs(p,t,f1,f2,f3,nu,alpha) or b=kstok3dp1brhs(p,t,f1,f2,f3,nu,alpha,vol,g) % [b,bp]=kstok2dp1brhs(p,t,f1,f2,nu,alpha) or [b,bp]=kstok2dp1brhs(p,t,f1,f2,nu,alpha,vol,g) %- ----------------------------------------------------------------------------

3D case