Sovann En

Alexis Lechervy

Frédéric Jurie

RPNet: an End-to-End Network for Relative Camera Pose Estimation

Keywords: relative pose estimation, pose estimation, posenet

This paper addresses the task of relative camera pose estimation from raw image pixels, by means of deep neural networks. The proposed RPNet network takes pairs of images as input and directly infers the relative poses, without the need of camera intrinsic/extrinsic. While state-of-the-art systems based on SIFT + RANSAC, are able to recover the translation vector only up to scale, RPNet is trained to produce the full translation vector, in an end-to-end way. Experimental results on the Cambridge Landmark data set show very promising results regarding the recovery of the full translation vector. They also show that RPNet produces more accurate and more stable results than traditional approaches, especially for hard images (repetitive textures, textureless images, etc.). To the best of our knowledge, RPNet is the first attempt to recover full translation vectors in relative pose estimation.

Introduction

In this paper, we are interested in relative camera pose estimation -a task consisting in accurately estimating the location and orientation of the camera with respect to another camera's reference system. Relative pose estimation is an essential task for many computer vision problems, such as Structure from Motion (SfM), Simultaneous Localisation And Mapping (SLAM), etc. Traditionally, this task can be accomplished by i) extracting sparse keypoints (ex. SIFT, SURF), ii) establishing 2D correspondences between keypoints and iii) estimating the essential matrix using 5-points or 8-point algorithms [START_REF] Nistér | An efficient solution to the five-point relative pose problem[END_REF]. RANSAC is very often used to reject outliers in a robust manner.

This technique, although it has been considered as the de facto standard for many years, presents two main drawbacks. First, the quality of the estimation depends heavily on the correspondence assignment. This is to say, too few correspondences (textureless objects) or too many noisy correspondences (repetitive texture or too much viewpoint change) can lead to surprisingly bad results. Second, the traditional method is able to estimate the translation vector only up to scale (directional vector).

In this paper, our objective is three folds: i) we propose a system producing more stable results ii) recovering the full translation vector iii) and we provide insights regarding relative pose estimation (i.e.. from absolute pose, from a pose regressor etc.).

As pointed out in [START_REF] Walch | Image-based localization using lstms for structured feature correlation[END_REF], CNN based methods are able to produce pretty good results in some cases where SIFT-based methods fail (i.e. texture less images). This is the reason why we opted for a global method based on CNN. Inspired by the success of PoseNet [START_REF] Kendall | Posenet: A convolutional network for realtime 6-dof camera relocalization[END_REF], we propose a modified Siamese PoseNet for relative camera pose estimation, dubbed as RPNet, with different ways to infer the relative pose. To the best of our knowledge, [START_REF] Melekhov | Relative camera pose estimation using convolutional neural networks[END_REF] is the only end-to-end system aiming at solving relative camera pose using deep learning approach. However, their system estimate the translation vector up to scale, while ours produces full translation vectors.

The rest of the paper is organized as follows: Section 2 presents the related work. Section 3 introduces the network architecture and the training methodology. Section 4 discusses the datasets and presents the experimental validation of the approach. Finally, Section 5 concludes the paper.

State of the Art

Local keypoint-based approaches. They address relative camera pose estimation using the epipolar geometry between 2D-2D correspondences of keypoints. Early attempts aimed at better engineering interest point detectors to focus on interesting image properties such as corners [START_REF] Harris | A combined corner and edge detector[END_REF], blobs in scale-space [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], regions [START_REF] Matas | Robust wide-baseline stereo from maximally stable extremal regions[END_REF], or speed [START_REF] Bay | Surf: Speeded up robust features[END_REF][START_REF] Tola | Daisy: An efficient dense descriptor applied to widebaseline stereo[END_REF][START_REF] Rublee | Orb: An efficient alternative to sift or surf[END_REF] etc. More recently, there is a growing interest to train interest point detectors together with the matching function [START_REF] Han | Matchnet: Unifying feature and metric learning for patch-based matching[END_REF][START_REF] Zagoruyko | Learning to compare image patches via convolutional neural networks[END_REF][START_REF] Tian | L2-net: Deep learning of discriminative patch descriptor in euclidean space[END_REF][START_REF] Detone | Superpoint: Self-supervised interest point detection and description[END_REF][START_REF] Trujillo | Using evolution to learn how to perform interest point detection[END_REF]. LIFT [START_REF] Yi | Lift: Learned invariant feature transform[END_REF] adopted the traditional pipeline combining a detector, an orientation estimator, and a descriptor, tied together with differentiable operations and learned end-to-end. [START_REF] Altwaijry | Learning to detect and match keypoints with deep architectures[END_REF] proposed a multitask network with different sub-branches to operate on varying input sizes. [START_REF] Detone | Superpoint: Self-supervised interest point detection and description[END_REF] proposed a bootstrapping strategy by first learning on simple synthetic data and increasing the training set with real images in a second time. End-to-End pose estimation. The first end-to-end neural network for camera pose estimation from single RGB images is PoseNet [START_REF] Kendall | Posenet: A convolutional network for realtime 6-dof camera relocalization[END_REF]. It is based on GoogLeNet with two output branches to regress translations and rotations. PoseNet followup includes: Baysian PoseNet [START_REF] Kendall | Modelling uncertainty in deep learning for camera relocalization[END_REF], Posenet-LSTM [START_REF] Walch | Image-based localization using lstms for structured feature correlation[END_REF] where LSTM is used to model the context of the images, Geometric-PoseNet where the loss is calculated using the re-projection error of the coordinates using the predicted pose and the ground truth [START_REF] Kendall | Geometric loss functions for camera pose regression with deep learning[END_REF]. Since all the 3D models used for comparisons are created using SIFT-based techniques, traditional approach seems more accurate. [START_REF] Walch | Image-based localization using lstms for structured feature correlation[END_REF] showed that the classical approaches completely fail with less textured datasets such as the proposed TMU-LSI dataset. [START_REF] Purkait | Spp-net: Deep absolute pose regression with synthetic views[END_REF] is an end-to-end system for pose regression taking sparse keypoint as inputs. Regarding relative pose estimation, [START_REF] Melekhov | Relative camera pose estimation using convolutional neural networks[END_REF] is the only system we are aware of. Their network is based on ResNet35 with FCs layers acting as pose regressor. Similar to the previous networks, the authors formulate the loss function as minimising the L2-distances between the ground truth and the estimated pose. Unfortunately, several aspects of their results (including Fig. 1. Illustration of the proposed system their label generation, experimental methodology and the baseline system) make comparisons difficult. Along side with pose regression problems, another promising works from [START_REF] Rocco | Convolutional neural network architecture for geometric matching[END_REF] showed that an end-to-end neural network can effectively be trained to regress to infer the homography between two images. Finally, two recent papers [START_REF] Yi | Learning to find good correspondences[END_REF][START_REF] Brachmann | Dsac-differentiable ransac for camera localization[END_REF] made useful contributions to the training of end-to-end systems for pose estimation. [START_REF] Yi | Learning to find good correspondences[END_REF] proposed a regressor network to produce essential matrix which can be then used to find the relative pose. However, their system is able to find the translation up to scale which is completely different from our objective. In [START_REF] Brachmann | Dsac-differentiable ransac for camera localization[END_REF], a differentiable RANSAC is proposed for outlier rejection and can be a plug-and-play component into an end-to-end system.

Relative pose inference with RPNet

Architecture. The architecture of the proposed RPNet, illustrated Fig. 1, is made of two building blocks: i) a Siamese Network with two branches regressing one pose per image, ii) a pose inference module for computing the relative pose between the cameras. We provide three variants of the pose inference module:

(1) a parameter-free module, (2) a parameter-free module with additional losses (same as PoseNet loss [START_REF] Kendall | Posenet: A convolutional network for realtime 6-dof camera relocalization[END_REF]) aiming at regressing the two camera poses as well as the relative pose, and (3) a relative pose regressor based on FC layers. The whole network is trained end-to-end for relative pose estimation. Inspired by PoseNet [START_REF] Kendall | Posenet: A convolutional network for realtime 6-dof camera relocalization[END_REF], the feature extraction network is based on the GoogLeNet architecture with 22 CNN layers and 6 inception modules. We only normalize the quaternion during test time. It outputs one pose per image.

For RPNet and RPNet + , the module for computing the relative pose between the cameras is straightforward and relies on simple geometry. Following the convention of OpenCV, the relative pose is calculated in the reference system of the 2nd camera. Let (R 1 , t 1 , R 2 , t 2 ) be the rotation matrices and translation vectors used to project a point X from world coordinate system to a fixed camera system (camera 1 & 2). Let (q 1 , q 2 ) be the corresponding quaternions of (R 1 , R 2 ). The relative pose is calculated as followed: where q * 1 is the conjugate of q 1 , and × denotes the multiplication in the quaternion domain. Both equations are differentiable. For RPNetFC, the pose inference module is a simple stacked fully connected layers with relu activation. To limit over-fitting, we modified the output of the Siamese network by reducing its output dimension from 2048 to 256. This results in almost 50% reduction of the number of parameters compared to PoseNet, RPNet and RPNet + network. The pose regressor network contains two FC layers (both with 128 dimensions). Losses. The loss function uses the Euclidean distance to compare predicted relative rotation q 1,2 and translation T1,2 with ground truth q1,2 and q 1,2 :

R 1,2 = q 2 × q * 1 and T 1,2 = R 2 (-R T 1 t 1 ) + t 2 (1)
loss = i (|| T i 1,2 -T i 1,2 || 2 + β * ||q i 1,2 -q i 1,2 || 2 ).
Quaternions are unit quaternions. The original PoseNet has a β term in front of quaternions to balance the loss values between the translation and rotation. To find the most suitable value of β, we cross-validated on our validation set. Please refer to our codes for different hyperparameter values on different subsets.

Experimentations

Experimental Setup

Dataset. Experimental validation is done on the Cambridge Landmark dataset 1 . Each image is associated with a ground-truth pose. We provide results on 4 of the 5 subsets (scenes). As discussed by several people, the 'street' scene raises several issues 2 . Pair generation. For each sequence of each scene, we randomly pair each image with eight different images of the same sequence. For a fair comparison with SURF, the pair generation is done by making sure that they overlap enough. We followed the train-test splits defined with the data set. Images are scaled so that the smallest dimension is 256 pixels, keeping its original aspect ratio. During training, we use 224*224 random crops and feed them into the network. During test time, we center crop the image. Baseline. The baseline is a traditional keypoint-based method (SURF). The focal length and the principle point are provided by the dataset. Other parameters are cross-validated on the validation set. For a fair comparison, we provide two scenarios for baselines: (1) the image are scaled to be 256*455 pixels, followed by a center-crop (224*224 pixels) to produce the same image pairs as tested with our networks and (2) the original images without down-sampling. We named Fig. 2. Translation and Rotation errors (median) of the different approaches these two scenarios as 'SURFSmall' and 'SURFFull'. All the camera parameters are adapted to the scaling and cropping we applied. Evaluation metric. We measured 3 different errors: i) translation errors, in meters ii) rotation errors, in degrees and iii) translation errors in degrees. We report the median for all the measurements.

Experimental results

Relative pose inference module. Fig. 2 compares the performance of the different systems and test scenarios. Based on these experimental results, RP-NetFC and RPNet + are the most efficient ways to recover the relative pose. On easy dataset (i.e.. KingsCollege and OldHospital), where there is no ambiguity textures, using pose regressor (RPNetFC) produces slightly better results than inferring the relative pose from the two images (PoseNet/RPNet/RPNet + ). On the contrary, on hard datasets (i.e.. ShopFacade and StMarysChurch), RPNetfamily outperforms RPNetFC. This behavior is also true for relative rotation and relative translation measured in degree. Globally, RPNetFC produces the best results followed by RPNet + , PoseNet and finally, RPNet. The differences of their results are between 0 and 8 degrees. Regarding technical aspect, RPNetFC is a lot easier to train than RPNet + /RPNet since it does not involve multiple hyper-parameters to balance the different losses. It also converges faster.

Comparison with traditional approaches. We will start by discussing the SURFSmall scenario first. In general, the error on both translation and rotation can be reduced between 5 to 70% using RPNet family, except on KingsCollege where the traditional approach slightly outperforms RPNet-based methods. We observed that the performance of the traditional approaches varies largely from one subset to another, while RPNet + /RPNetFC are more stable. In addition, the traditional approach requires camera information for each image in order to correctly estimate the pose. In contrast, RPNet-based does not require any specific information at all. Using the original image size (SURFFull) significantly boost the performance of the traditional approach. However, RPNetFC still enjoy a significant gain in performance on OldHospital and ShopFacade, while performing slightly worse than SURFFull on KingsCollege and StMarysChurch. The difference in performance between SURFFull and RPNetFC is even more significant when the images contain large view point changes (see Fig. 3). Full translation vector. One of our objectives is to provide a system able to estimate the full translation vector. On average, we observed that the median error ranges between 2 to 4 meters, using RPNetFC. Fig. 4 gives an idea of ground truth translations w.r.t. reference axes (xyz). For instance, on KingsCollege, the values of X-axis can range from -29m to 30m with an STD of 7 meters. Interestingly, our network has a translation error of only meters.

Conclusions

This paper proposed a novel architecture for estimating full relative poses using an end-to-end trained neural network. The network is based on a Siamese architecture, which was experimented with different ways to infer the relative poses. In addition, to produce competitive or better results over the traditional SURF-based approaches, our system is able to produce an accurate full translation vector. We hope this paper will provide more insight and motivate other researchers to focus on global end-to-end system for relative pose regression problems.

Fig. 3 .

 3 Fig. 3. Accumulative hist. of errors in rotation (1st row, d), translation (2nd row, m).

Fig. 4 .

 4 Fig. 4. Min/Max/Mean/STD relative translations (ground truth), w.r.t. XYZ axis (m).

Table 1 .

 1 Number of training and testing pairs for Cambridge Landmark dataset. SE stands for spatial extent, measured in meter.

	Scene	Train Test SE	Scence	Train Test SE
	Kings College 9.1k 2.4k 140x40 Shop Facade	1.6k 0.6k 35x25
	Old Hospital 6.5k 1.2k 50x40 St Marys Church 11k 4.1k 80x60

http://mi.eng.cam.ac.uk/projects/relocalisation

https://github.com/alexgkendall/caffe-posenet/issues/2

Acknowledgements. This work was partly funded by the FrenchUK MCM ITP program and by the ANR-16-CE23-0006 program.


F Chinesta
email: chinesta@ec-nantes.fr

A Poitou

S Chatel

J A Garcia

L L Gascon

ON THE EQUIVALENT IN-PLANE PERMEABILITY

Keywords: RTM, Permeability, Separated representations, Model reduction

Simulation of RTM processes is usually performed in 2D. Thus, only a mesh of the middle plane is needed with the associated degrees of freedom savings with respect to fully 3D modeling. However such modeling needs the definition of an equivalent in-plane permeability representing the ignored dimension (the thickness). The definition of such permeability is not a trivial task because each ply in the thickness direction can be anisotropic, being the principal anisotropy direction different from one ply to the neighbor plies. In this work we propose a novel fully 3D modeling whose computational cost is equivalent to a 2D solution. It allows addressing properly the equivalent in-plane permeability issue.

INTRODUCTION

In general the simulation of RTM processes assumes a 2D flow model. The most usual model results of combining the Darcy's law and the flow incompressibility:

0 p          v K v
that results in the second order BVP:

  0 p     K
The main issue in defining this model concerns the definition of the permeability tensor K . Different techniques exist, but in that follows we are assuming that an averaged permeability has been determined for each type of reinforcement architecture. First we are considering a laminate composed of many layers with the same principal directions, whose permeability tensor writes 0 0 0 0

               

K

Now, we consider the flow in a plate mould, being the permeability principal directions parallel to the mould walls. In this case the permeability tensor becomes diagonal. The simplest scenario concerns    and 0   . In this case the flow is unidirectional. When the ply is placed after performing a certain in-plane rotation to the principal permeability directions, the  component of the permeability tensor is no more zero. In that case when we impose a constant injection pressure, being the pressure in the opposite mould boundary zero, the resulting pressure field is depicted in Fig. 1. We can conclude that if the principal directions of the permeability of each layer are aligned along the coordinate axes, even if the permeability evolves in the thickness direction, the flow remains unidirectional. The situation is radically different if the permeability also evolves in the x-direction. When we assumes an evolution of   , x z  as depicted in figure 5, the pressure field that results, illustrated in figure 6, implies the velocity field depicted in figure 7, which exhibit noticeable deviations with respect to an unidirectional flow. These simple examples have proven that in general laminates, involving many layers with different reinforcement orientations the situation could be quite complex and the assumption of a simple average of the permeability of the different layers could be a crude approximation that deserves a deeper analysis. For this purpose a fully 3D simulation could be the best alternative, but an accurate representation of the thickness direction could become too expensive from the computational view point. An appealing alternative consists in applying a separated representation of the pressure field involved in the porous media flow modelling. This kind of approximation allows computing 3D solutions with a computational cost characteristic of 2D simulations, that is, of standard RTM simulations. In the next section we summarized the main elements of such technique, deeply described in many of our former works [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representations[END_REF], and that is here particularized to the treatment of the Darcy's model defined in plate moulds. Its extension to shell geometries is quite direct.

PGD IN PLATE DOMAINS

In what follows we are illustrating the construction of the Proper Generalized Decomposition of a model defined in a plate domain

I    with 2    and   0, I H  :   0 p     K (1)
We consider that the laminate is composed of P different anisotropic plies each one characterized by a well defined permeability tensor   , i x y K -it is assumed constant in the ply thickness-. Moreover, without a loss of generality, we assume the same thickness for the different layers of the laminate, that we denotes by h. Thus, we can define a characteristic function representing the position of each layer:

1 1 ( ) 0 i i i z z z z otherwise         , 1, , i P  L (2)
where  

1 i z i h    . Now, the laminated
permeability can be given in the following separated form:

      1 , , i P i i i x y z z       K K x (3)
      1 , i N i i i p z X Z z      x x (5)
In what follows we are illustrating the construction of one such decomposition. For this purpose we assume that at iteration n N  the solution is already known:

      1 , i n n i i i p z X Z z      x x (6)
and that at the present iteration we look for the solution enrichment:

    1 , , ( ) ( )
n n p z p z R S z     x x x (7)
The test function involved in the weak form is searched under the form:

  * * * , ( ) ( ) ( ) ( ) p z R S z R S z     x x x (8)
By introducing Eqs. ( 7) and ( 8) into (4) it results:

°°°°°* * * * * * * * n R S R S R S d dS dS dS R R R dz dz dz R S R S d dS dS R R dz dz                                                                                                           K Q ( 9 )
where ° denotes the plane component of the gradient

operator °  , T x y       and n Q denotes the flux at iteration n: °        1 i i i n n i i i X Z z dZ z X dz                   x Q K x ( 1
    1 n R X   x x and     1 n S z Z z   .
We are illustrating each one of the just referred steps:

1. Computing   R x from   S z : When   S z is known the test function reduces to:   * * , ( ) ( ) p z R S z   x x (11)
and the weak form (9) reduces to:

°°°* * * * n R S R S d dS dS R R dz dz R S d dS R dz                                                             K Q (12)
Now, as all the functions involving the coordinate z are known, they could be integrated in

  0, I H  . Thus, if
we consider:

T         k K k k ( 1 3 ) with xx x y yx yy k k k k        k , xz yz k k        k and zz k   , then we can define: 2 0 0 2 0 0 z H z H z z x z H z H T z z dS S dz Sdz dz dS dS Sdz dz dz dz                                   k K k k (14) and °    0 0 1 0 0 z H z H i i i N z z x n i z H z H i i T i i z z dZ SZ dz S dz dz X X dZ dS dS Z dz dz dz dz dz                                                       k x Q x k k (15) that allows writing equation (12) into the form °°°* * * * x x n R R d R R R d R                                         K Q (16)
that defines an elliptic 2D problem defined in the middle plane of the plate.

Computing

 

S z from   R x : When   R x is known the test function writes:   * * , ( ) ( ) p z R S z   x x (17)
and the weak form (9) reduces to:

°°°* * * * n R S R S d dS dS R R dz dz R S d dS R dz                                                              K Q (18)
Now, as all the functions involving the in-plane coordinates   , x y  x are known, they could be integrated in  . Thus, using the previous notation, we can define:

°  °    °  °    2 z R R d R Rd R Rd R d                                     k K k k (19) and °°  °°1 i i i i N z n i i i i R X d R Xd Z dZ dz R X d RX d                                                              k Q k k ( 20
                                                K Q (21)
that defines a one-dimensional BVP.

CONCLUSIONS

This paper comes back to a recurrent issue in the numerical modelling of RTM flows, the one related to the pertinence of using 2D flow models making use of an averaged permeability of the different layers involved in the laminate. In complex situations significant deviations could be found, these deviations are being studied at present, and could justify the use of a fully 3D modelling. However, 3D simulations are reputed expensive from the point of view of the computational resources required for addressing complex scenarios. The use of separated representations as the ones involved in the proper generalized decompositions -PGD-could be an appealing alternative for addressing 3D models with a cost characteristic of 2D simulations. The application of the PGD on the fully 3D simulation of flows encountered in RTM processes constitutes a work in progress.

Figure 1 :

 1 Figure 1: Pressure distribution

Figure 2 :

 2 Figure 2: Evolution of the permeability in the thickness

Figure 3 :Figure 4 :

 34 Figure 3: Pressure distribution

Figure 5 :

 5 Figure 5: Evolution of the permeability

  z is searched under the separated form:

  0 ) Now, as the enrichment process is non-linear we propose to search the couple of functions   R x and   S z by applying an alternating direction fixed point algorithm. Thus, assuming   R x known, we compute   S z , and then we update   R x . The process continues until reaching convergence. The converged solutions allow defining the next term in the finite sums decomposition:


Ll Gascón

J A García
email: jugarcia@mcm.upv.es

F Chinesta

E Ruiz

F Trochu

ANALYSIS OF THE SATURATION IN LIQUID COMPOSITE MOLDING PROCESSES USING AN ESSENTIALLY NON-OSCILLATORY (ENO) TECHNIQUE

Keywords: Liquid Composite Molding, Saturation, Void, Essentially Non-Oscillatory Techniques

In the manufacturing of composite parts by Liquid Composite Molding process (LCM), complete saturation of the fibrous reinforcement is key. Incomplete saturation leads to voids within the fibers which cause failure of the final product. Thus, understanding of the formation of voids is necessary for proper molding. In order to analyze the formation of voids during the resin impregnation process, a one-dimensional solution based on two-phase flow through a porous medium, is proposed. This model leads to the introduction of relative permeability as a function of saturation and a modified equation for the saturation as a non-linear advection-diffusion equation with viscous and capillary phenomena, which depends on a number of factors. A detailed analysis is performed to assess the relative significance of the various input parameters on the saturation profiles. In order to numerically solve the modified saturation equation for the LCM process a high order essentially non-oscillatory (ENO) technique is proposed. The implemented algorithm allows a numerical optimization of the injected flow rate, which minimizes the micro/macroscopic void formation during mold filling. Some numerical results are presented and compared with the results taken from literature in order to validate the proposed mathematical model and the numerical scheme.

INTRODUCTION

Equations that describe the LCM filling process with void formation are based on a two phase flow model and lead to a coupled system of a nonlinear advection-diffusion equation for saturation and an elliptic equation for pressure and velocity [START_REF] García | An efficient solver of the saturation equation in liquid composite molding processes[END_REF]. In general, the saturation equation is a non-linear advection-diffusion equation which includes the capillary pressure effect and it reduces to a purely advection transport equation when capillary effects are neglected. The hyperbolic nature of the saturation equation and its strong coupling through relative permeability represent a challenging numerical issue. In this paper, a very accurate numerical approach is proposed to solve this complex flow behavior. In previous works [START_REF] Breard | Numerical simulation of void formation in LCM[END_REF][START_REF] Ruiz | Optimization of injection flow rate to minimize micro/macro-voids formation in resin transfer molded composites[END_REF] different experiments were carried out to investigate the process of void formation. These experimental works shown that macrovoids tend to form during injection at low flow rates, due to capillary dominant effects, whereas high injection rates lead to microvoids formation. Based on these observations, a simplified mathematical model has been proposed in order to model the LCM filling process with void formation. In this model, the diffusivity coefficient in the saturation equation has been replaced by a term which depends on the velocity. The permeability is assumed to be a function of saturation, and then the continuity equation that governs the pressure distribution, includes a source term which depends on the saturation. Essential to the optimum process design in LCM is the numerical simulation of the modified saturation equation. Many numerical methods to solve this type of equations suffer from serious nonphysical oscillations, excessive numerical dispersion or a combination of both. The technique here used for solving the advection-diffusion equation which governs the evolution of the degree of saturation of porous media is based on an essentially non-oscillatory fixed mesh strategy. For the ENO schemes, interpolation polynomials of one order less than the order of accuracy required in the solution are computed and these polynomials are a good approximation to the values of the numerical flux function at the cell walls. The key idea in the rth-order ENO schemes is to use the "smoothest" stencil among r possible candidates to approximate the fluxes at cell boundaries to high-order accuracy and at the same time to avoid spurious oscillations near shocks. However, they also have certain drawbacks. One problem is with the freely adaptive stencil choice, which could change.

GOVERNING EQUATIONS

The mathematical formulation of the saturation in LCM takes into account the interaction between resin and air as it occurs in a two phase flow. Combining equations that describe mass conservation and Darcy's laws for resin and air phases as described in [START_REF] García | An efficient solver of the saturation equation in liquid composite molding processes[END_REF], the resulting equation for the saturation in its most general form gives S + ( ) ( )

cf vf S D S S t (1)
where

( ) , ( )
( ) R c cf A R A S P f S D S f S S S S S (2)
Here ( ) cf D S is the nonlinear diffusivity coefficient due to capillary pressure P c , defined as c A R P P P ; v is the total velocity, S is the degree of saturation of the reinforcement by the liquid resin, ( ) ( ) / j Rj j S K S K is the phase mobility, with ( ) Rj K S the relative permeability of the phase j, j the viscosity of phase j and K the permeability tensor. Replacing the total velocity v by R A v v and simplifying, Eqn. 1 can be rewritten as follows

S + ( ) 1 ( ) R R R P v f S f S S S t S ( 3
)
Assuming that the diffusivity coefficient depends on the resin flow velocity, as identified in experimental observations, leads to the following simplified equation,

S + ( ) M R m R R v f S v S t v (4)
The term on the right side has been proposed in [1 M m represent the dispersive coefficients of the macro and microscopic voids, respectively. Since the main interest is to simulate the flow of the resin phase, we can derive the governing equations for this problem combining, for the resin phase, Darcy's law (the subscript R for v and p has been omitted), ( )

R s a t K S K v p (5) with mass conservation S v t (6)
To derive a closed model, Eqn. 4 for the saturation has been considered. The simulation of the filling process involves the following operations at each time step: 1. Calculate the pressure distribution by applying a standard finite element discretization to Equation S ( )

R sat K S p K t ( 7
)
where the relative permeability and the term on the right side depend on the saturation degree. 2. Calculate the velocity field from Darcy's law for the resin. 3. Update the saturation distribution by integrating Eqn. 4 using a fourth-order ENO technique (Algorithm ENO-4). The boundary conditions are given by: the pressure gradient in the normal direction to the mold walls vanishes, the pressure or the flow rate is specified on the inflow boundary and the pressure is zero in the empty part of mold.

NUMERICAL SCHEME FOR THE SATURATION

In this study we describe a fourth-order essentially non-oscillatory scheme for the numerical saturation solution. The reconstruction algorithm is based on an adaptive selection of stencil for each cell so as to avoid spurious oscillations near discontinuities while achieving high order of accuracy away from them. For the saturation equation S + ( )

v f S D S t (8)
we define the flux as 1 ( )

F vf S D S (9)
and consider a uniformly spaced grid where each cell and then taking a derivative. Following the technique described in [START_REF] Fedkiw | The Penultimate Scheme for Systems of Conservation Laws: Finite Difference ENO with Marquina's Flux Splitting[END_REF], the calculation of H(x) is based on polynomial interpolation. The main ingredient of the ENO method is the adaptive choice of stencil. For our case, it begins with a starting point to the left or right of the current cell by means of upwinding determined by the sign of the velocity; as the order of the undivided differences is increased, the undivided differences themselves determine the stencil: the smaller undivided difference is chosen from two possible choices at each stage, ensuring a smoothest fit. Then the following fourth-order algorithm based on the ENO-Roe numerical flux reconstruction can be described.

ENO-4 Algorithm for the saturation:

Step 1. Compute

1 1 1 2 1 2 3 1 1 4 2 1 1 2 3 1 1 ( ) , , , 2 2 1 1 2
, 3 3 6 24

j j j j j j j j j j j j j j j j j v v v DF F DF F F h D F F F F D F F F F F h h Step 2.
If D D j k k r r s v k j k j F F r k r k F F s r s r F then 1, else s F t s t s 1
Step3. Compute

2 2 1 1 2 3 2 3 3 3 4 + 2 1 3 1 4 6 2( ) 2 k r s j t F DF h D F j k h D F j r h D F j s j s t s
In the numerical simulation of the modified saturation equation two terms contribute to smooth the flow front: one is related to the source term, and the other purely numerical term is introduced by the dicretization scheme. The last effect can be reduced by using a fourth-order numerical scheme. Fig. 1 shows numerical results using the ENO-4 Algorithm for the saturation of the onedimensional test described in [START_REF] García | An efficient solver of the saturation equation in liquid composite molding processes[END_REF]. A mold of length 1 m has been considered. The saturated permeability K sat and the resin viscosity are set to 10 8 m 2 and 0.1 Pa.s, respectively. For the numerical simulation, we take Fig. 1 Saturation profiles for a 1D RTM filling at constant flow rate [START_REF] García | An efficient solver of the saturation equation in liquid composite molding processes[END_REF]

  injection rate of 0.001 m/s. The domain is assumed initially empty, except the first element that represents the injection nozzle that is assumed full-filled. Two different models for Eqn. 8 has been illustrated: ( ) f S S , at the left, Vel=-0.001 m/s. A=1e-010. a=0.001. Void%= 12.Vel=-0.001 m/s. A=1e-010. a=0.001. Void%= 2.3918 Lenght, m (N° Iterations=3601) Saturation, S MFA= 37.2658%

ACKNOWLEDGMENTS

This research work is supported by a grant from the Ministerio de Educación y Ciencia (MEC), project DPI2007-66723-C02-01.


Ilaria Fragal À

Bozhidar Velichkov

SERRIN-TYPE THEOREMS ON TRIANGLES

Keywords: overdetermined problems, triangles, torsional rigidity, principal frequency, pcapacity. 2010 Mathematics Subject Classification: 35N25, 35J57, 49Q10

.

Introduction

Starting from the celebrated symmetry result by Serrin [START_REF] Serrin | A symmetry problem in potential theory[END_REF], a huge literature is devoted to the study of overdetermined boundary value problems of various kind, both for interior and exterior problems, and for a large variety of differential operators and overdetermined conditions. With no attempt of completeness, see [START_REF] Bianchini | An overdetermined problem for the anisotropic capacity[END_REF][START_REF] Brandolini | Serrin-type overdetermined problems: an alternative proof[END_REF][START_REF] Brock | A symmetry result for an overdetermined elliptic problem using continuous rearrangement and domain derivative[END_REF][START_REF] Brock | Some new symmetry results for elliptic problems on the sphere and in Euclidean space[END_REF][START_REF] Crasta | On the Dirichlet and Serrin problems for the inhomogeneous infinity Laplacian in convex domains: regularity and geometric results[END_REF][START_REF] Fragalà | Symmetry results for overdetermined problems on convex domains via Brunn-Minkowski inequalities[END_REF][START_REF] Fragalà | Partially overdetermined elliptic boundary value problems[END_REF][START_REF] Fragalà | Overdetermined problems with possibly degenerate ellipticity, a geometric approach[END_REF][START_REF] Garofalo | A symmetry result related to some overdetermined boundary value problems[END_REF][START_REF] Kawohl | Overdetermined problems and the p-Laplacian[END_REF][START_REF] Reichel | Radial symmetry by moving planes for semilinear elliptic BVPs on annuli and other non-convex domains, Elliptic and parabolic problems (Pont-à-Mousson[END_REF][START_REF] Reichel | Radial symmetry for elliptic boundary-value problems on exterior domains[END_REF][START_REF] Shahgholian | Diversifications of Serrin's and related symmetry problems[END_REF][START_REF] Silvestre | Overdetermined problems for fully nonlinear elliptic equations[END_REF][START_REF] Vogel | Symmetry and regularity for general regions having a solution to certain overdetermined boundary value problems[END_REF] and references therein.

In this setting we feel a bit surprising that no symmetry result is available, to the best of our knowledge, for overdetermined problems set on polygons. This seems to be an interesting topic especially from the point of view of Calculus of Variations. In fact, polygons which support solutions to certain overdetermined problems turn out to be stationary, in the sense of shape derivation, for some some long-standing open problems in shape optimization; among these, let us mention the minimization of the principal frequency among polygons with a given number of sides under an area constraint, for which the regular gon is conjectured to be optimal [START_REF] Henrot | Extremum problems for eigenvalues of elliptic operators[END_REF]Section 3.3], or the minimization of the p-capacity among convex planar sets under a perimeter constraint, for which the expected solution is a degenerate polygon, i.e. a line segment [START_REF] Bucur | Optimal convex shapes for concave functionals[END_REF]Remark 2.4].

In this paper we give a first contribution to the study of overdetermined problems on polygons, by considering the case of triangles. As a model case, we start from the torsion problem as in Serrin's paper, namely we consider the Dirichlet problem

-∆u = 1 in T, u = 0 on ∂T, (1.1)
where we fix the space dimension to be n = 2, and T lies in the class T of (nondegenerate) triangles. We overdetermine problem (1.1) with one of the following setting of boundary conditions of integral type for |∇u| on ∂T :

A i+1 A i i 2 -|x -A i | |∇u(x)| 2 dx = 0 ∀i ∈ {1, 2, 3} , (1.2)
A i+1 A i i 2 -|x -A i | |∇u(x)| 2 dx = κ i 2 f (θ i ) -f (θ i+1 ) ∀i ∈ {1, 2, 3} , (1.3)
where A i are the vertices of T , i is the length of the side [A i A i+1 ], κ is a positive constant, and

f (θ) := cot θ + 1 sin θ . (1.4)
Conditions (1.2) and (1.3) are very natural candidates to yield symmetry. Indeed, as we will show in Section 2, they correspond precisely to the stationarity conditions for the torsional rigidity respectively under a volume and a perimeter constraint, when one rotates a side of T around its midpoint. Since the equilateral triangle is known to be the unique maximizer of torsional rigidity among triangles of given volume or perimeter, it is reasonable to expect it may be also the unique "critical" triangle.

These arguments can be repeated if the torsion functional is replaced by the principal frequency, namely we may consider in place of (1.1) the Dirichlet problem

-∆u = λ 1 (T )u in T, u = 0 on ∂T, T u 2 dx = 1, (1.5)
where λ 1 (T ) is the first Dirichlet eigenvalue of the laplacian.

For both problems (1.1) and (1.5), we prove the following Serrin-type result:

Theorem 1.1 (triangular symmetry, interior case). Let T ∈ T , and let u be the unique solution to (1.1) or to (1.5). Then:

-u satisfies the overdetermined boundary condition (1.2) if and only if T is equilateral.

-u satisfies the overdetermined boundary condition (1.3) for some κ ≥ 0 if and only if T is equilateral.

Our approach is based on a reflection argument, and is reminiscent of Serrin's proof in the use of Hopf's boundary point principle. In particular, similarly as shown by Reichel in [START_REF] Reichel | Radial symmetry for elliptic boundary-value problems on exterior domains[END_REF] for Serrin's result, our symmetry statement can be extended to the case of exterior problems. However, this extension is a partial one, in the sense that it only concerns the overdetermined boundary condition (1.2). To be more precise, for a fixed exponent p ∈ (1, 2), denote by ∆ p the p-laplace operator, and consider the p-capacitary problem

∆ p u = 0 in R 2 \ T, u = 1 on ∂T, lim |x|→∞ u(x) = 0. (1.6)
We have: 3) holds only on one side of the triangle, then the angles adjacent to such a side must be equal, so that the triangle is isosceles with the selected side as a basis. Thus, in order to obtain that T is equilateral, it is enough to impose the overdetermined condition (1.2) (respectively, (1.3)) on two sides of the triangle. This refinement can be interpreted as a symmetry result for a partially overdetermined boundary value problem, in the spirit of [START_REF] Fragalà | Partially overdetermined elliptic boundary value problems[END_REF]. As a second case of study, we overdetermine problems (1.1), or (1.5), or (1.6), with a different setting of boundary conditions, which looks as well quite natural, as they correspond to the stationarity conditions for the associated energy functional, respectively under a volume and a perimeter constraint, when one makes a side of T move in a parallel way to itself (cf. Section 2). They read

A i+1 A i |∇u(x)| 2 dx = κ i ∀i ∈ {1, 2, 3} , (1.7)
A i+1 A i |∇u(x)| 2 dx = κ f (θ i ) + f (θ i+1 ) ∀i ∈ {1, 2, 3} , (1.8)
where κ is a positive constant, f (θ) is given by (1.4), and as above we intend that the term |∇u(x)| 2 is changed into |∇u(x)| p when dealing with problem (1.6).

In contrast to conditions (1.2)-(1.3), conditions (1.7)-(1.8) do not yield symmetry. The difference is in fact even more drastic, since the last setting of conditions turns out to be always satisfied. We have indeed: Theorem 1.6 (triangular equidistribution). Let T ∈ T , let u be the unique solution to (1.1), (1.5), or (1.6), and let F be respectively the torsional rigidity, the first Dirichlet eigenvalue of the Laplacian, and the p-capacity. Then

-u satisfies (1.7), with κ = |α| 2
F (T ) Vol(T ) -u satisfies (1.8), with κ = |α| F (T ) Per(T ) ,
where α is the homogeneity degree of F under domain dilation (divided by (p -1) in case of problem (1.6)).

Only in the case of the first Laplacian Dirichlet eigenvalue and just for condition (1.7), Theorem 1.6 was proved by Christianson in the recent paper [START_REF] Christianson | Equidistribution of Neumann data mass on triangles[END_REF], by using a completely different approach. Our proof seems to be more direct, and has the advantage to work the same way in all the cases covered by the statement, and possibly even in more general situations, such as for instance the cases of p-torsion function and first Dirichlet eigenvalue of the p-Laplacian. (Avoiding the technicalities related to more sophisticated versions of maximum principles is the reason why we preferred to present problems (1.1) and (1.5) in their linear version).

To extend our results to quadrilaterals or polygons with an arbitrary number of sides seems to require some different idea and remains by now a challenging open problem.

Stationarity conditions

In this section we consider scale invariant energies defined on the class K 2 of convex bodies with nonempty interior in R 2 , and we show that the equalities (1.2)-(1.3) and (1.7)-(1.8) correspond respectively to the stationarity conditions for such energies, when the perturbed domain is a triangle and the perturbation consists in rotating one side around its midpoint and in moving it in a parallel way to itself.

To be more precise, we fix the following setting:

Definition 2.1 (energy functionals). Let Ω ∈ K 2 . Let F be either the torsional rigidity, or the first eigenvalue of the Laplacian, or the p-capacity (1 < p < 2), defined respectively as

τ (Ω) := -inf u∈H 1 0 (Ω) Ω |∇u| 2 -2 u dx , λ 1 (Ω) := inf u∈H 1 0 (Ω)\{0} Ω |∇u| 2 dx Ω |u| 2 dx , Cap p (Ω) := inf R 2 |∇u| p dx : u ∈ Ẇ 1,p (R 2 ) , u = 1 on Ω ,
where Ẇ 1,p (R 2 ) denotes the closure of space of the space of smooth compactly supported functions with respect to the norm R 2 |∇u| p dx. We denote by α the homogeneity degree of F under dilations, so that

α =      4 if F = τ -2 if F = λ 1 2 -p if F = Cap p .
Moreover, we denote by u Ω respectively either the torsion function, or the first Dirichlet eigenfunction of the Laplacian, of the p-capacitary potential of Ω, namely the solution to the Euler-Lagrange equation of F given in (1.1), (1.5), (1.6).

Next we introduce the two family of deformations we are going to consider for triangles. Definition 2.2 (rotation around the midpoint). Given T ∈ T with vertices A 1 , A 2 and A 3 , and t ∈ R small enough, we denote by {T t } the triangle with vertices A t 1 , A t 2 , A t 3 , obtained by keeping fixed A 3 =: A t 3 and replacing A 1 and A 2 by A t 1 and A t 2 in such a way that (see Figure 1): • the line determined by A t 1 and A t 2 is a rotation on an angle t of the line determined by A 1 and A 2 around the midpoint M of the side [A 1 , A 2 ]; • A t 1 and A t 2 lie respectively on the lines determined by A 3 , A 1 and A 3 , A 2 ; • the direction of the rotation is determined by the convention that, for t > 0, A t 2 lies on the segment [A 2 , A 3 ], while for t < 0, A t 1 lies on the segment

[A 1 , A 3 ]. T A 1 A 2 A 3 = A t 3 A t 2 A t 1 M 0 > t T A 1 A 2 A 3 = A t 3 A t 2 A t 1 0 < t Figure 1.
The rotating around the midpoint (on the left) and the parallel movement of a side (on the right). We are now ready to identify stationary triangles:

Definition 2.3 (parallel movement). Given T ∈ T with vertices

Proposition 2.4. Let E Vol (T ) and E Per (T ) be defined by The above statement is obtained as a direct consequence of Lemmas 2.6 and Lemma 2.7 below. In turn, to obtain such results, we need a more general first derivation formula, which is valid for suitable perturbations of convex domains, as stated in Lemma 2.5 below. Let us recall that, for any Ω ∈ K 2 , the gradient of the function u Ω in Definition 2.1 is well defined on ∂Ω and belongs to L p (∂Ω) [START_REF] Dahlberg | Estimates of harmonic measure[END_REF][START_REF] Lewis | Regularity and free boundary regularity for the p Laplacian in Lipschitz and C 1 domains[END_REF]. Here and throughout the remaining of the paper, we mean that the exponent p equals 2 when dealing with the interior problems (1.1) and (1.5).

E Vol (Ω) := F (Ω) 1 α Vol(Ω) 1/2 , E Per (Ω) := F (Ω)

Lemma 2.5 (Shape derivatives with respect to generic perturbations).

Let Ω ∈ K 2 , with unit outer normal n Ω , and let

Ω t = Φ t (Ω), where t ∈ [0, T ) → Φ t ∈ W 1,∞ (R 2 ) is differentiable at t = 0, with Φ 0 (x) = x and d dt t=0 Φ t (x) = V (x). Let F and u Ω be as in Definition 2.1. Then d dt t=0 F (Ω t ) = (sign α) (p -1) ∂Ω V • n Ω |∇u Ω | p dH 1 . (2.2)
Proof. When F = τ , the statement can be deduced by applying the duality approach developed in [START_REF] Bouchitté | Shape derivatives for minima of integral functionals[END_REF] for shape derivative of minima of integral functionals (cf. in particular [2, Theorems 3.7 and Example 3.9], and see also [START_REF] Bouchitté | A variational method for second order shape derivatives[END_REF]Section 7.2]). When F = λ 1 , we refer to [19, Section 5.9.3]. When F = Cap p , one can adapt the duality approach in [START_REF] Bouchitté | Shape derivatives for minima of integral functionals[END_REF], by showing that [2, Theorems 3.3 and 3.7] continue to hold for the exterior problem of p-capacity and give directly (2.2). However, for the sake of completeness, we give here a formal derivation of (2.2), by sketching a proof along the more classical line adopted in [START_REF] Fragalà | On an isoperimetric inequality for capacity conjectured by Pólya and Szegő[END_REF] to compute the shape derivative ot the 2-capacity of convex bodies in R 3 .

Let u t be the p-capacitary potential on Ω t , that is

-div(|∇u t | p-2 ∇u t ) = 0 in R 2 \ Ω t , u t = 1 on ∂Ω t .
By arguing in a similar way as in the Appendix of [START_REF] Fragalà | On an isoperimetric inequality for capacity conjectured by Pólya and Szegő[END_REF], we see that the map R t → u t ∈ W 1,p loc (R 2 \ Ω) is differentiable in t = 0. Denoting by u the derivative of u t at t = 0 and formally differentiating the equation for u t , we get that u is a solution of

-div |∇u| p-2 ∇u + (p -2)|∇u| p-4 (∇u • ∇u )∇u = 0 in R 2 \ Ω, u = -V • ∇u on ∂Ω. (2.3) Let now ψ ∈ C ∞ c (R 2
) be a smooth function such that ψ ≡ 1 in a neighborhood of Ω. Testing (2.3) with u -ψ and integrating by parts we get

Ω c |∇u| p-2 ∇u • ∇(u -ψ) dx = -(p -2) Ω c |∇u| p-4 (∇u • ∇u )(∇u • ∇(u -ψ)) dx = -(p -2) Ω c |∇u| p-2 (∇u • ∇u ) dx + (p -2) Ω c |∇u| p-4 (∇u • ∇u )(∇u • ∇ψ) dx, which gives (p -1) Ω c |∇u| p-2 ∇u • ∇u dx = Ω c |∇u| p-2 ∇u • ∇ψ + (p -2)|∇u| p-4 (∇u • ∇u )(∇u • ∇ψ) dx.
(2.4) On the other hand, denoting by n t the unit outer normal to Ω t , we have

Cap p (Ω t ) = Ω c t |∇u t | p dx = - ∂Ωt |∇u t | p-2 ∂u t ∂n t ψ dH 1 = Ω c t |∇u t | p-2 ∇u t • ∇ψ dx.
Taking the derivative in t = 0, exlpoiting (2.4) and integrating by parts, we obtain

d dt t=0 Cap p (Ω t ) = Ω c |∇u| p-2 ∇u • ∇ψ + (p -2)|∇u| p-4 (∇u • ∇u )(∇u • ∇ψ) dx = (p -1) Ω c |∇u| p-2 ∇u • ∇u dx = -(p -1) ∂Ω |∇u| p-2 (∇u • n Ω )u dH 1 = (p -1) ∂Ω |∇u| p-2 (∇u • n Ω )(V • ∇u) dH 1 ,
which is precisely (2.2), since ∇u = |∇u| n Ω on ∂Ω. Actually let us recall from [START_REF] Lewis | Regularity and free boundary regularity for the p Laplacian in Lipschitz and C 1 domains[END_REF] that the gradient ∇u Ω (y) has (non-tangential) limits as y → x ∈ ∂Ω for H 1 -a.e. x ∈ ∂Ω (see also [START_REF] Colesanti | The Hadamard variational formula and the Minkowski problem for p-capacity[END_REF]Lemma 2.13]), and that, denoting such limits by ∇u Ω (x), we have |∇u Ω | ∈ L p (∂Ω, dH 1 ).

Lemma 2.6 (Shape derivatives with respect to side rotations). Let F be as in Definition 2.1.

Let T ∈ T and let {T t } be as in Definition 2.2. Then

d dt t=0 F (T t ) = (sign α) (p -1) M A 1 |x -M | |∇u T | p dH 1 (x) - A 2 M |x -M | |∇u T | p dH 1 (x) = (sign α) (p -1) A 2 A 1 1 2 -|x -A 1 | |∇u T | p dH 1 (x), (2.5)
d dt t=0 Vol (T t ) = 0, (2.6)
d dt t=0 Per (T t ) = 1 2 (f (θ 1 ) -f (θ 2 )), where f (θ) = cot θ + 1 sin θ . (2.7)
Proof. If Φ t : R 2 → R 2 is a one-parameter family of diffeomorphisms such that Φ t (T ) = T t , then by an elementary geometic argument we see that for every

x ∈ [A 1 , A 2 ] such that |x -M | < 1 2 min 1, 1 + sin t cos θ 1 sin(θ 1 -t) , 1 - sin t cos θ 2 sin(θ 2 + t) ,
we have Vol (T t ) = 1 .

Φ t (x) • n T (x) =        |x -M | sin t sin θ 1 sin(θ 1 -t) 1 + sin t cos θ 1 sin(θ 1 -t) -1 if x ∈ [A 1 , M ], -|x -M | sin t sin θ 2 sin(θ 2 + t) 1 - sin t cos θ 2 sin(θ 2 + t) -1 if x ∈ [M, A 2 ]. Then V (x) := d dt t=0 Φ t (x) satisfies V (x) • n T (x) = |x -M | if x ∈ [A 1 , M ], -|x -M | if x ∈ [M, A 2 ].
(2.9)

d dt t=0 Per (T t ) = f (θ 1 ) + f (θ 2 ).
(2.10)

Proof. Notice that the perturbation described in Definition 2.3 corresponds to the diffeomorphism Φ t (x) = h+t h x, where without loss of generality we suppose A 3 = 0 and we set h to be the distance from

A 3 to the side [A 1 , A 2 ]. Thus V (x) = x/h and V (x) • n T (x) ≡ 1 on [A 1 , A 2 ]
, so we get (2.8), while (2.9) and (2.10) follow directly by the equalities

Per (T t ) = Per (T ) + t f (θ 1 ) + f (θ 2 )
and Vol (T t ) = Vol (T ) + t 1 + o(t).

3. Proof of Theorems 1.1 and 1.2

We will prove Theorems 1.1 and 1.2 by a reflection argument. We first give the geometric construction, which is the same for the cases of interior or exterior domains, and then we will treat them separately. Geometric construction. Suppose by contradiction that T has two different inner angles θ 1 < θ 2 . Let M be the midpoint of the side [A 1 , A 2 ] and let L M be the line passing through M and orthogonal to the side [A 1 , A 2 ]. Without loss of generality we can suppose that M = 0, L M is the y-axis {x = 0}, A 1 = (-1 /2, 0), A 2 = ( 1 /2, 0) and A 3 = (x 3 , y 3 ) with x 3 > 0 and y 3 > 0 (see Figure 2). Let N be the intersection of L M with the side [A 1 , A 3 ]. We denote by Ω int the interior of the triangle with vertices M, N and A 2 and by Ω ext the unbounded domain {x > 0} \ T . For any function u : R 2 → R we consider the reflection ũ : R 2 → R of u defined by ũ(x, y) = u(-x, y).

A 1 A 2 A 3 M N L M θ 1 θ 2 Ω int Ω ext

Figure 2. Construction of the reflected domain

Proof of Theorem 1.1. We carry out the argument for the torsion function, namely the solution u T to (1.1), being the case of the first eigenfunction of the Dirichlet Laplacian completely analogous. By elliptic regularity, u T is C 1,α up to the boundary of T (except at the vertices), and by the strong maximum principle we have that u T > 0 in T (see for instance [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]). We now consider the reflection ũT of u T . We notice that ∆(u T -ũT ) = 0 in Ω int , while on the boundary ∂Ω int we have:

ũT = u T on [N, M ] ∪ [M, A 2 ] and 0 = ũT < u T on [N, A 2 ].
Thus, u T -ũT > 0 in Ω int and by Hopf's boundary point lemma

|∇u T | > |∇ũ T | on [M, A 2 ]. (3.1)
Multiplying this inequality by |x -M | and integrating on the segment [M, A 2 ], we get

A 2 M |x -M | |∇u T | 2 dH 1 (x) > A 2 M |x -M | |∇ũ T | 2 dH 1 (x) = M A 1 |x -M | |∇u T | 2 dH 1 (x) , (3.2)
which implies

A 2 A 1 1 2 -|x -A 1 | |∇u T | 2 dH 1 (x) < 0 . (3.3)
This is clearly in contradiction with (1.2), and also with (1.3) since θ 1 < θ 2 and the function θ → f (θ) is monotone decreasing on (0, π).

Proof of Theorem 1.2. Let now u T : R 2 \ T → R be the p-capacitary potential of T . By a reflection argument with respect to each side, we get that u T is C 1,α up to the boundary of T (except at the vertices). By the strong maximum principle for p-harmonic functions, we have the strict inequality u T < 1 on the open set Ω ext = R 2 \T (see [START_REF] Manfredi | p-harmonic functions in the plane[END_REF]). We notice that ∆ p ũT = 0 = ∆ p u T in Ω ext , while on the boundary ∂Ω ext we have: By the strong comparison principle for p-harmonic functions (see [START_REF] Manfredi | p-harmonic functions in the plane[END_REF]) we get that ũT < u T in Ω ext . As a consequence, by Hopf's principle for p-harmonic functions (see for instance Let T ∈ T . We observe that, if {T t } is the one-parameter family of triangles given by Definition 2.3, for every t the triangle T t is homothetic to T . Then, since the functionals E Vol (T ) and E Per (T ) from The derivatives appearing in the above equations can be explicitly computed by using formulas (2.8)-(2.9)-(2.10) in Lemma 2.6. This way we see that the two equalities (4.1) correspond exactly to (1.7) and (1.8), the value of κ being given as in the statement of Theorem 1.6. Since T is an arbitrary triangle, our proof is achieved.

Remark 1 . 5 .

 15 Problem (1.6) overdetermined by (1.3) (with |∇u| 2 replaced by |∇u| p ) remains interestingly open. This discrepancy with respect to the overdetermined problems (1.1)-(1.3) and (1.5)-(1.3) covered by Theorems 1.1 should not surprise the attentive reader. Indeed, while condition (1.2) represents a "discrete" analogous to the stationarity condition under volume constraint which for smooth domains reads |∇u| = c, condition (1.3) should be assimilated to the stationarity condition under perimeter constraint which for smooth domains reads |∇u| = c(H), being c a function of the boundary mean curvature H. This latter overdetermined condition, treated by Serrin in [26,Theorem 3], can be successfully handled by the moving planes method only under a specific monotonicity assumption on c, and the favourable sign of the monotonicity changes when passing from interior to exterior problems.

  ũT = u T on L M ∪ [M, A 2 ] and ũT < u T = 1 on [N, A 3 ] ∪ [A 2 , A 3 ].

[23, Chapter 5 , 2 M 2 M 2 A 1 1 2 4 .

 522214 Section 5.5]) we obtain|∇ũ T | > |∇u T | on [M, A 2 ].(3.4)Then, similarly as in the case of interior problems, integrating this inequality we getM A 1 |x -M | |∇u T | p dH 1 (x) = A |x -M | |∇ũ T | p dH 1 (x) > A |x -M | |∇u T | p dH 1 (x),(3.5)which impliesA -|x -A 1 | |∇u T | p dH 1 (x) Proof of Theorem 1.6

  (2.1) are invariant under dilations, it holds d dt t=0 E Vol (T t ) = 0 and d dt t=0 E Per (T t ) = 0 . (4.1)

  A 1 , A 2 and A 3 , and t ∈ R small enough, we denote by {T t } the triangle with vertices A t 1 , A t 2 , A t 3 , obtained by keeping fixed A 3 =: A t 3 and replacing A 1 and A 2 by A t 1 and A t 2 in such a way that (see Figure ... ): • the line determined by A t 1 and A t 2 is parallel to the side [A 1 , A 2 ] and at distance |t| from [A 1 , A 2 ]; • A t 1 and A t 2 lie respectively on the lines determined by A 3 , A 1 and A 3 , A 2 ; • the direction of the movement is determined by the convention that, for t > 0, A t 2 and A t 1 do not lie to the segment [A 2 , A 3 ] and to the segment [A 1 , A 3 ] respectively.

  Vol (T ) and E Per (T ) with respect to the perturbations {T t } as in Definition 2.2, if and only if u T satisfies respectively (1.2) and (1.3). (ii) A triangle T ∈ T is stationary for E Vol (T ) and E Per (T ) with respect to the perturbations {T t } as in Definition 2.3, if and only if u T satisfies respectively (1.7) and (1.8).

	1
	α Per(Ω)	,	(2.1)

with F as in Definition 2.1.

(i) A triangle T ∈ T is stationary for E

  Lemma 2.7 (Shape derivatives with respect to parallel movements). Let F be as in Definition 2.1. Let T ∈ T and let {T t } be as in Definition 2.3. Then

	d dt t=0	F (T t ) = (sign α) (p -1)	A 1 A 2	|∇u| p dH 1 .	(2.8)
			d
			dt t=0
	Now, applying Lemma 2.5 to this family Φ t , we get (2.5).
	On the other hand, through elementary geometric arguments, we see that
	Per (T t ) = Per (T ) -1 +	1 sin θ 1 2 sin(θ 1 -t)	+	1 sin θ 2 2 sin(θ 2 + t)	+	1 sin t 2 sin(θ 1 -t)	-	1 sin t 2 sin(θ 2 + t)
	Vol (T t ) = Vol (T ) + o(t),
	which immediately give (2.6) and (2.7).

Acknowledgments. The authors have been supported by the Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM). The second author has been partially supported by ANR through the projects GeoSpec (LabEx PERSYVAL-Lab, ANR-11-LABX-0025-01) and CoMeDiC (ANR-15-CE40-0006).

We are grateful to the referee for pointing out the refinement discussed in Remark 1.4.