A DMLPG Refinement Technique for 2D and 3D Potential Problems

Annamaria Mazzia,Giorgio Piniand Flavio Sartoretto

A DMLPG Refinement Technique for 2D and 3D Potential Problems

Annamaria Mazzia1,Giorgio Pini1and Flavio Sartoretto2

Meshless Local Petrov Galerkin(MLPG)methods are pure meshless techniques for solving Partial Differential Equations(PDE).MLPG techniques are nowadays used for solving a huge number of complex,real-life problems.While MLPG aims to approximate the solution of a given differential problem,its “dual”Direct MLPG(DMLPG)technique relies upon approximating linear functionals.Assume adaptive methods are to be implemented.When using a mesh-based method,inserting and/or deleting a node implies complex adjustment of connections.Meshless methods are more apt to implement adaptivity,since they does not require such adjustments.Nevertheless,ad-hoc insertion and/or deletion algorithms must be devised,in order to attain a good accuracy.In this paper we introduce a fresh refinement technique for DMLPG methods.Nodes are inserted in a discretization cloud where the local variation in the solution is supposed to be“large”.The variation is estimated using the(local)Total Variation(TV).DMLPG allows to directly estimate the partial derivatives,in order to compute the TV.MLPG must rely upon approximating the derivatives of the shape functions,hence MLPG refinement results to be more involved than its DMLPG counterpart.We show that our DMLPG refinement procedure allows one to efficiently solve a given diffusion problem whose solution undergoes large variations on a small portion of the domain.The accuracy afforded by a fine uniform cloud can be attained by using far less non-uniformly arranged nodes.

Meshless Methods;MLPG;Direct MLPG;Refinement;Diffusion.

1 Introduction

Meshless methods nowadays supply a good alternative or integrative technique to Finite Element(FE)methods.As an example,meshless methods apply as an eli-gible choice for fracture analysis,see Yang,Budarapu,Mahapatra,Bordas,Zi,and Rabczuk(2015).

Meshless methods are particularly attractive when meshing is a time-consuming and cumbersome task[Gerace,Erhart,Kassab,and Divo(2013)],as it happens in adaptive methods.

Kriging was proposed as a method to identify a suitable discretization cloud[Chen and Liew(2011)].This is a quite involved procedure.A key element for adaptivity,is a refinement procedure,which allows one to refine a discretization only where the solution u of the given differential problem undergoes large variations.In this paper,we introduce a refinement procedure,based on approximating the variation in u by estimating the(local)Total Variation(TV)of the numerical solution˜u.

Meshless Petrov Galerkin(MLPG)methods[Atluri(2004)]are truly meshless,and widely used.Our previous paper[Mazzia,Pini,and Sartoretto(2014)]deals with devising a refinement procedure for MLPG,that is based on the Moving Least Squares(MLS)method proposed in Lancaster and Salkauskas(1981).

Direct MLPG(DMLPG)[Mirzaei and Schaback(2013)],exploits the Generalized Moving Least Squares method(GMLS)[Mirzaei,Schaback,and Dehghan(2012)].While MLS aims at approximating a multivariate function,GMLS approximates linear functionals.The solution of a linear Partial Differential Equation(PDE),formulated in strong or weak form,can be computed by approximating via GMLS the associated functionals.DMLPG exploits GMLS to approximate weak forms of PDE problems.DMLPG can be seen as a “dual”method to MLPG.

To identify the domain regions where any refinement is required,we adopted the(local)Total Variation(TV)as an indicator of the local variation in u.Estimating the TV entails computing the partial derivatives.DMLPG enables us to approximate the partial derivatives directly,while other meshless methods,MLPG included,must exploit expensive approximations of the partial derivatives of u.

To analyze the numerical properties of a refinement method,one must consider sound test problems and ad-hoc,simple domains.In the sequel,we numerically analyze the accuracy and efficiency achieved with our DMLPG-based refinement procedure when solving Poisson problems on either[0,1]2or[0,1]3.

This paper is organized as follows.Section 2 recalls basic facts about the approximating schemes exploited by MLPG and DMLPG.Section 3 sketches the fundamental concepts underlying MLPG and DMLPG techniques.Section 4 identifies the trial and test spaces,which are a key ingredient of our numerical methods.Section 5 details our refinement strategy.Section 6 describes our distinguished test problems.Section 7 describes and discusses our numerical results.Section 8 summarizes our conclusions.

2 Approximation schemes

The Moving Least Squares(MLS)method[Lancaster and Salkauskas(1981)]has been proposed for approximating a function,u(),inside a region Ω⊂Rdgiven a number,N,of its values,In the context of MLPG methods,in order to approximate a given d-dimensional problem,MLS is exploited for generating a set of d-dimensional trial functions on the grounds of a suitable “weight”function[Mazzia,Pini,and Sartoretto(2008);Mazzia and Sartoretto(2010)].The ensuing trial functions are called the“shape”functions of MLS.

The Generalized Moving Least Squares(GMLS)method[Mirzaei,Schaback,and Dehghan(2012)],is a technique which generalizes the MLS formulation[Levin(1998)],aiming to approximate continuous linear functionals in the dual of(Ω),for any k≥ 0.GMLS aims to approximate a linear functional,λ(u),based on a given set of N linear functionals λi(u).One obtains an approximation

Each coefficient φi(λ ),called a GMLS shape function,must be linear in λ .As an example,λ(u)can be a partial derivative of u,whilecan be a given set of u values on given nodesin this case can be interpreted as a generalization of Finite Difference methods,based on allvalues,instead of a given stencil.

For a given continuous linear functional λ(u),the polynomial-based GMLS computesan approximationwhere p∗minimizes a weighted least-squares error functional[Mazzia,Pini,and Sartoretto(2012)].

For more details on the algorithm,see e.g.Levin(1998);Mirzaei and Schaback(2013).

3 MLPG and DMLPG techniques

Let us consider the linear Poisson equation on the domain Ω

where f is a given source function,being any point in Ω.Dirichlet and Neumann boundary conditions are imposed on the domain boundary∂Ω

In order to compute an approximationof the solution,a finite set of trial functionsis chosen.

Let us assume that the residual of eq.(2)is multiplied by a suitable test function τ.The divergence theorem is applied,thus obtaining the weak formulation for(2)

for any τ in a suitable functional spaceMany MLPG methods have been proposed[Atluri(2004);Fries and Matthies(2004)],each of which can be identified by an appropriate choice of trial and test functions.

In order to approximate the solution of our weak formulation,a set of discretization nodes must be given.Let N be the total number of nodes.

A set of trial functions,ξi,is usually adopted,each of which is“centered”on nodeThe support ofis a ball centered atwhose radius is ri.Actually,for each i we setoutside Ω,which is equivalent to consideringinstead ofThe support of each test functionis a ball(or a cuboid)centered atwhose radius(or half side length)is ρi.We assume τi=0 outside Ω,thus consideringWe assume that τi=0 is ona typical setting in many MLPG schemes[Mirzaei and Schaback(2013);Mazzia and Sartoretto(2010)].In principle the test functions need not necessarily be centered on theFor the sake of simplicity,we use as support centers thenodes,for both the trial and the test functions.

A set of Local Weak Forms(LWF)is obtained by writing eq.(4)for each test function

The DMLPG technique is implemented by applying GMLS to the weak problem(5).The solution of our Poisson problem is approximated by using a polynomial space.

All the details of our implementation of the DMLPG technique for diffusion problems are given in Mazzia,Pini,and Sartoretto(2012).

4 Finite dimensional spaces

On the grounds of our previous 2D and 3D results[Mazzia,Pini,and Sartoretto(2008);Mazzia and Sartoretto(2010);Mazzia,Pini,and Sartoretto(2014)],we exploit appropriate weight functions w(t)[Atluri and Zhu(2002);Belytschko,Krongauz,Organ,Fleming,and Krysl(1996);Lancaster and Salkauskas(1981);Mazzia,Pini,and Sartoretto(2008)],in order to identify suitable trial and test spaces.

4.1 Generating functions

Let us assume that w(t)is a given,differentiable,compact supported,generator function;when t＞1,w(t)=0 holds.

where w(t)is a Gaussian generator according to Lu,Belytschko,and Gu(1994),i.e.

Note that σ is a parameter controlling the function shape.Throughout this paper we assume that σ=1.

Our test functions τiare Tensor Product Functions(TPF).We give the de finitions for 3D problems;2D problems are treated by disregarding the z values.With each node,we associate the TPF test function

The function f(t)is in turn a polynomial generator according to Liu(2009)

4.2 Identifying the radiuses

A crucial step in obtaining accurate(D)MLPG schemes involves identifying the support trial basis radiuses ri,and the test basis ones ρi,i=1,...,N[Mazzia and Sartoretto(2010)].

For each node we arrange in ascending order its distances from the closest nc≪N neighbors.The ncinteger value was identified by suitable numerical experiments(see below).

Let us assume a given irregular discretization I.

In the case of 2D domains,for each nodein I we compute thenc=7nodes closest toWe arrange their distances in ascending order

(a)If dk=d,k=1,...,4,we can assume that the node distribution around the i-th node is“practically uniform”.We set ri=βd,ρi=αd,α,β are parameters to be tuned.See the sequel for details.

(b)If not,we set

The boundary nodes are treated in the same way as the internal nodes.

The parameters α and β were identified by means of numerical experiments:They fall in not too large ranges[Mazzia and Sartoretto(2010)].Typical intervals are

Throughout this paper we assume that α=1.This setting proved suitable for all the test problems shown in the sequel.

As for 3D problems,it should be noted that there are six faces of a cube centered on any given pointWe perform step(a)on six points,i.e. “local uniformity”corresponds to

Alternatively,we perform step(b)by assuming that nc=9.

As we saw before,the parameter setting α=1 was adopted.

5 Discretizations

5.1 Uniform grids

Let us assume for the sake of simplicity that our domain is either the[0,1]2square,for 2D problems,or the[0,1]3cube for 3D problems.

When dealing with 2D problems,let us start with the uniform grid on[0,1]2identified by intersecting the sides of the square and those three x-parallel and three y-parallel lines that divide the x-and y-side evenly into nx=ny=4 parts.Let us call this uniform discretization U1,the “level”ℓ=1 uniform discretization.The interval spacing is h1=1/4,N=5×5=25 nodes are identified.Each finer uniform discretization level is obtained by halving each sub-interval.We thus obtain N=25,81,289,1089,4225 nodes,corresponding to levels=1,2,3,4,5.

The 3D uniform discretization grids on[0,1]3are obtained in a similar manner,by setting nx=ny=nz=4.On the initial level ℓ=1,we thus have h1=1/4,N=5×5×5=125 nodes are enrolled.Each finer level is obtained by halving each sub-interval.We obtain N=125,729,4913,35937,274625 nodes,corresponding to levels N=1,2,3,4,5.

5.2 Refinement strategy

The efficacy for DMLPG computations of the irregular discretizations built using our refinement strategy has to be tested.

For any irregular discretization,there are two measures worth considering.They are the fill distance hI,Ω,and the separation distance,qI,i.e.Fasshauer(2007)

5.2.1 2D problems

The first discretization level=1 is set to I1=U1,the N=25 node uniform grid.To re fine a given discretization Ikat level=k,we consider on each nodethe“local”Total Variation(TV)of the solution u,de fined as

where uxand uyare the partial derivatives of the solution.The advantages of this“variation measure”are clearly explained in Strang,(2007),pag.415:The variational methods which use the TV norm“allow for discontinuities but disfavor oscillations”.

To approximate the TV on nodeby taking our numerical solutioninto account,we set

|Ωi|being the area of the i-th DMLPG integration subdomain.One can prove that,due to our previous assumptions on the domains of the trial and test functions,Ωiis the support of the i-th test function,i.e.the square centered atthe radius of which is ρi.

Figure 1:Node insertion 2D strategy,when a“locally-uniform”discretization is detected around node

Let

If we assume a given threshold parameter γ,we re fine our discretization around each given nodeif and only if

Our refinement procedure around nodeinvolves adding one node in the middle of each line joiningwith thenodes closest to it.The valuemust be guessed on the grounds of geometrical considerations.We set=8 for 2D problems,so that all“diagonal”nodes are added on an uniform grid(see Figure 1).

Note that when a node being inserted overlaps an old one,its insertion is skipped.Once all the nodes in a given discretization Ikhave been processed,we say that a new “discretization level”Ik+1is ready.

Our refinement procedure is repeated until a maximum number Nmaxof nodes has been reached,or a maximum number of“discretization levels”has been computed.The maximum number of nodes depends on the storage space available.

Note that,under the previous assumptions,uniform refinements are obtained when γ=0 and an initial uniform distribution has been adopted,On the other hand,uniform refinements can be easily computed by evenly dividing intervals,whereas performing our refinement procedure demands some more assessments.That is why,although the nodes are the same,the CPU time spent on computing an uniform refinementby setting γ=0 is considerably longer than the time needed to compute the nodes directly inThis applies to the 3D domain,[0,1]3too.

5.2.2 3D problems

The first discretization level=1 is set to I1=U1,the N=125 node uniform grid on[0,1]3.

When dealing with 3D problems,the straightforward generalization of our 2D“local”TV applies,where subdomain volumes are involved instead of areas.

Our refinement strategy is updated accordingly by setting=24.When a(local)uniform discretization is detected,all the “diagonal”nodes in a cube.As in 2D domains,the setting γ=0 produces uniform clouds.

5.3 Implementation issues

We implemented our algorithm in FORTRAN 77 codes,compiled via XLF v9.1.They were run on a machine operating under Ubuntu with an Intel Core i7-2600K,3.40 GHz(quad core)processor,and 2x4GB,1333 MHz,RAM.

No complicated data structures were needed to deal with our discretizations because they are mere clouds of points.

Our Fortran codes use a CSR sparse matrix storage representation to deal with the large number of nodes needed for volume discretizations in non-structural,(e.g.fl ow,heat)problems.

To arrange in ascending order the distances of each node from its neighbors,we used the Fortran subroutine SORT2 available in the NAPACK library.

6 Test problems

To check our adaptive strategy,we assign the forcing function f and compute the boundary conditions in eq.(2),so that its“test”solution is a function u undergoing large variations on a small portion of the domain.

6.1 2D problems

First,we consider the classical Gaussian function,centered on a given point P0=(x0,y0),i.e.

Figure 2:Contour regions for the solution of the test problem P2(T).

The parameter c is a large positive value that generates a high “hump”around P0.In the sequel,we set c=200 unless stated otherwise.

Let us assume that we numerically solve the Poisson problem(2)in Ω=[0,1]2,having set the Dirichlet boundary conditions such that its solution is the function(12).The setting P0=(1/2,1/2),the centroid of our domain,corresponds to the 2D problem called P2(GC)in the sequel,where“GC”stands for“Gaussiancentroid”.

Any adaptive procedure is likely to be effective when finer discretizations adopt a large number of discretization nodes near the point P0where a large variation in u occurs.On the other hand,“far”away from P0the u values are small,and u does not display large variations,so the nodes can be distributed quite coarsely with no appreciable loss of accuracy.

As a further test problem we consider,as in Kee,Liu,Zhang,and Lu(2008)

This function displays a “hill”rising from the bottom left of[0,1]2.Figure 2 shows the contour levels of the surface.

Let us assume that we numerically solve the Poisson problem(2)in Ω=[0,1]2,having set the Dirichlet boundary conditions such that its solution is the func-tion(13).

The ensuing differential problem is labeled test problem P2(T),where“T”stands for the“arcTan-based”solution.

6.2 3D problems

Simply extending our 2D test problems,we consider Gaussian functions centered at a given point P0=(x0,y0,z0)

When dealing with 3D problems,we set c=200 unless stated otherwise.

Let us assume that we aim to solve the Poisson problem on[0,1]3with Dirichlet boundary conditions such that the function(14)is its solution.We exploited many settings,corresponding to Gaussian “centers”on suitable domain points,either inside or on the boundary of our domain.

·The setting P0=(1/2,1/2,1/2),is called test problem P3(GC).The label“GC”stands for“Gaussian-Centroid”.

·When we set P0=(1/2,1/2,0),we are speaking of test problem P3(GXY).The label“GXY”reminds us that the Gaussian is centered on X=Y=1/2.

·The setting P0=(1/2,0,0),is called test problem P3(GX).The label“GX”stands for“Gaussian centered on X=1/2”(and Y=Z=0).

·The setting P0=(1/2,1/2,0),c=500,is called test problem P3(GXY5).The “GXY5”is a mnemonic for the GXY case,where “c=500 was set”.

We also consider the function

By applying the Dirichlet boundary conditions to the 3D Poisson equation,the ensuing test problem is called P3(T).The label“T”stands for“arcTan”.

7 Numerical results

N being the number of discretization nodes,we define the numerical error

Table 1:Number of nodes, fill distanceand separation distancecomputed for uniform and re fined 2D clouds,when solving problem P2(GC).The refined clouds were obtained by setting γ=10-3.

Table 1:Number of nodes, fill distanceand separation distancecomputed for uniform and re fined 2D clouds,when solving problem P2(GC).The refined clouds were obtained by setting γ=10-3.

Uniform γ=10-3ℓN hUℓ,[0,1]2 qUℓN hRℓ,[0,1]2 qRℓ1 25 1.77E-001 1.25E-001 25 1.77E-001 1.25E-001 2 81 8.84E-002 6.25E-002 81 8.84E-002 6.25E-002 3 289 4.42E-002 3.13E-002 289 4.42E-002 3.13E-002 4 1089 2.21E-002 1.56E-002 497 4.42E-002 1.56E-002 5 4225 1.10E-002 7.81E-003 1790 4.42E-002 7.81E-003

where uiis the exact value on nodewhileu˜iis the corresponding approximated value.

When solving 2D problems,we exploited cubic polynomials in the GMLS approximation.When solving 3D problems,we exploited the quadratic polynomial approximation space.

We are mainly interested in comparing the accuracy achieved by means of our refinement procedure as opposed to uniform discretizations.

7.1 2D Numerical Results

We performed several 2D numerical experiments,but to summarize here our 2D results we only report here on two test problems.Our procedure actually proves most effective on 3D problems,which demand a large number of nodes in order to be solved accurately.

Let us assume that β=3.5 for 2D problems.

Based on our previous experiments with MLPG refinements[Mazzia,Pini,and Sartoretto(2014)],and extensive numerical experiments,we set the refinement parameter γ in the range of 0 ≤ γ≤ 10-2.

7.1.1 Discretizations

For reasons of storage and efficiency,our finest uniform cloud is the one with N=4,225 nodes obtained at level=5.This is the finest uniform cloud that we used in all our 2D numerical experiments.

Table 1 shows the number of nodes,and the fill and separation distance for some discretization levelsgenerated when solving the test problem P2(GC).We consider the uniform cloudsand the corresponding refined onesobtained by setting γ=10-3.When using other γ values,and dealing with the other test prob-lems,the distances recorded come between the “uniform”and the γ=10-3types of behavior.

Figure 3:Cloud nodes adopted by setting γ=10-3at the refinement level ℓ=5,when Problem P2(T)is solved.

Figure 4:Errors recorded when our 2D test problems are solved.The errors obtained using uniform clouds coincide with the γ=0 refinement levels.Frame(a)refers to Problem P2(GC).Frame(b)applies to Problem P2(T),where the curve for γ=0 overlaps those corresponding to γ=1e-4,1e-3.

As we can see from the columns labeled“Uniform”in Table 1,at each new level the fill and separation distances are halved when an uniform refinement is performed,as expected.The refinement with γ=10-3produces the same separation distances,while the fill distanceremains unchanged at levels=3,4,5.Recalling the de finitions(11)one can argue that our refinement procedure generates“ fine sub-regions”the distances of which are larger than their internal separation,whilethe uniform refinement does not.

Table 2:The cloud level needed for attaining the HUA is given for some values of the refinement parameter γ.The corresponding number of nodes N,and the error e obtained are also shown.

We could show that the refined(non-uniform)clouds cluster around the Gaussian center P0,but we prefer to show the more interesting behavior seen when dealing with Problem P2(T).Figure 3 shows the nodes identified by our refinement procedure.The refinement level=5 is considered.By comparing Figure 3 with Figure 2,one can see that the nodes coalesce well around the “hill”corresponding to the region of greatest variation in the solution.

Recalling our considerations on the cloud distance measures,we argue that our refinement technique effectively clusters the nodes around the regions of high variation in the solution.

7.1.2 Convergence history

Figure 4 plots the behavior of the errors vs the number of nodes in the discretizations.The refinement parameter γ=0,1e-4,1e-3,1e-2 values were set.Looking at Figure 4,one can see that the error is almost always lower when the number of refined nodes is larger.Frame(b)in Figure 4 shows that the error increases slightly when we go formU1toU2,and this applies to non-uniform refinements too.

Figure 4 shows that,by setting γ=0,we can replicate the uniform clouds.Larger and larger values γ values produce non-uniform clouds that enable the same accuracy to be achieved with a smaller number of nodes.

Table 2 reports the number of refined nodes needed to to obtain the same accuracy order as with our finest,uniform discretization level,which contains N=4225 nodes.In the sequel this is called the “Highest Uniform Accuracy”(HUA).

It is easy to see by comparing Figure 4 with Table 2 that our refinement procedure enables the HUA to be reached using a considerably smaller number of nodes thanin our finest uniform cloud.

Table 3:Fill distances and separation distances computed for uniform and refined 3D clouds when solving problem P3(GX),γ=5×10-4.

As an example,row 4 in Table 2 reports that when dealing with problem P2(GC),by setting γ=1.0e-2,only N=1172 nodes in the refined cloud enable the HUA to be reached,which corresponds to N=4225 nodes in the uniform discretization.

We note that,as one can infer from Table 2,when dealing with the non-Gaussian solution(Problem P2(T)),more nodes are required to reach the HUA.This comes as no surprise:Gaussian solutions undergo large variations in the vicinity of the Gaussian center,and small variations elsewhere,while the“arcTan-based”test solution undergoes large variations over a much broader portion of the domain.

7.2 3D Numerical Results

By means of extensive numerical experiments we tuned 3.0≤β≤4.0,using the value that afforded the best accuracy for each test problem.

Our analysis and numerical experiments showed that 0≤γ≤5×10-3.

7.2.1 Discretizations

Table 3 shows the fill and separation distance for our 3D clouds generated when solving test problem P3(GX).The level=1 corresponds to an uniform N=125 node cloud,as mentioned earlier.We consider the uniform clouds and our refined clouds obtained by setting γ=10-3.The finest uniform cloud contains N=35937 nodes.Table 3 shows that the refined clouds have the same separation distances as the uniform grids.The fill distances in the uniform grids halve at each level,as one can guess.On the other hand,the fill distances in the non-uniform,refined clouds remain practically unchanged for levels ℓ=2,3,4.These results confirm that,as in 2D problems,our 3D refinement procedure generates fine sub-clouds containing more nodes than in other regions where a coarser discretization is produced.

Figure 5 sketches the cloud obtained with our refinement procedure at level=5,when Problem P3(GX)is dealt with.The refinement parameter value is γ=5×10-4.As expected,many nodes are clustered around the Gaussian center P0=(1/2,0,0).

Figure 5:Non-uniformly refined clouds at levels ℓ=5,6,when dealing with Problem P3(GX).The Gaussian center is P0=(1/2,0,0)The refinement parameter value was set to γ=5.0×10-4.

Figure 6:Errors recorded when our 3D test problems are solved,using some γ values.We recall that,by setting γ=0,one generates uniform clouds.Frame(a)refers to Problem P3(GC),and Frame(b)to Problem P3(GXY).

Note that the cloud obtained at level ℓ=6 produces exactly the same Figure.The refinement procedure adds points close to the “blob”shown,so in Figure 5 they overlap with the nodes at level ℓ=5.

Table 4:For each problem and refinement parameter value,we report the level ℓand the number of nodes enrolled for reaching the HUA(see the error values e).The CPU seconds Tγspent on solving the Problem and the“speedup”with respect to the time spent when using the finest uniform discretization TUare also shown.

These considerations,supplemented with our previous analysis on Table 3,confirm that our 3D refinement procedure enables an effective clustering of the nodes around the regions where the greatest variations in the solution occur.

7.2.2 Convergence history

Figure 7:Same as in the previous Figure.Frame(c)refers to Problem P3(GX),and Frame(d)to Problem P3(GXY5).

Figure 8:Errors recorded when Problem P3(T)is attacked.Note that the curves overlap almost completely.

Figure 9:Number of nodes in the HUA 3D discretization level vs the refinement parameter value γ.

Figures 6,7,and 8 report the errors raised when various refinement parameter values are set.Note that the error almost always decreases when the refining nodes are added to the discretization cloud.By inspecting these Figures,one can confirm that,like for 2D problems,by setting γ=0 one can replicate the uniform clouds,whereas setting γ＞ 0 produces non-uniform clouds.

We aim to attain the highest accuracy afforded by the finest uniform mesh(HUA),counting N=35,937 nodes.The HUA can clearly be achieved with a far smaller number of nodes,however,by exploiting our refinement procedure.

Table 4 summarizes our results.As shown in the Table,the larger the γ,the smallerthe number of nodes required to attain the HUA.By setting γ=5e-3 the HUA is attained on Gaussian-based problems P3(GC),P3(GXY5),P3(GXY),P3(GX),using less than 5200 nodes.The best case concerns Problem P3(GX),where N=1315≪35937 nodes in the refined non-uniform cloud suffice to arrive at the HUA.Now let us consider the curve pertaining to Problem P3(T),whose solution is“arcTan-based”.When γ=5e-3,Table 4 shows that N=30,165 nodes are required to reach the HUA.Such a large number is due to the large region where the solution of problem P3(T)undergoes high variations.

Table 5:FEM errors when the level ℓ=5 uniform mesh of Type(a)or Type(b)is used.To reach the HUA,ℓ=5 levels were generated.The finest uniform meshes adopt the same N=274,625 nodes(though their connections differ).

Figure 9 plots the number of discretization nodes needed to attain the HUA in relation to the refinement parameter γ.We recall that γ=0 corresponds to uniform discretizations.By inspecting this Figure,it easy to see that our refinement is very effective for Gaussian solutions,but less for the arcTan-based.

Let us compare our results with those obtained using a standard numerical method.Figure 11 shows the errors raised using linear tetrahedral Finite Elements vs the number of nodes.The plotted values were obtained using two types of 3D uniform mesh,called “Type(a)”and “Type(b)”in the sequel.Level ℓ=2 of our Type(a)and Type(b)uniform meshes are plotted in Figure 10.Our two meshes differ in the orientation of some sides of the triangles forming the base of some of the tetrahedral elements discretizing the inside of the domain cube.Note the different behavior of the errors in Frames(a)and(b)of Figure 11.One can see that using the same uniform distributed nodes,but choosing the elements differently,gives rise to a very different error behavior.This is a known FEM effect.The difference is exacerbated when the Gaussian centroid is set on a boundary side of the domain cube(see the curve for Problem P3(GX),X=1/2,Y=Z=0 in Figure 11).

These changes are not found when different clouds obtained by our refinement procedure are used.This suggests that our meshless procedure is less prone to node distribution in the discretization clouds than FEM is to mesh connections.

Figure 10:Sketch of the boundary triangles which form the bases of the tetrahedral elements inside one of our 3D FEM meshes.Frame(a)shows a mesh of Type(a),Frame(b)refers to Type(b).

Table 5 shows the FEM errors.To reach the HUA by FEM,the level ℓ=5 mesh,containing N=274,625 nodes,must be exploited.This number of mesh nodes is far larger than in our refined clouds;we needed to use less than N=35,937 nodes,which is the number of nodes in the ℓ=4 uniform grid(see Table 4).

7.2.3 Efficiency

Figure 11:Errors obtained when our 3D test problems are solved via linear FEM using uniform,tetrahedral meshes.Frame(a)shows the results obtained by Type(a)meshes.Frame(b)refers to Type(b)meshes.

Let TUbe the CPU seconds spent on computing the solution using our finest uniform cloud.We recall that performing our refinement procedure with the setting γ=0 produces the same nodes as in the uniform clouds.But the CPU time taken up by the former procedure is longer than when using“pre-computed”uniformclouds.When reporting CPU times,we must therefore distinguish between situations when“pre-computed”uniform discretizations are used as opposed to computing refinements with γ=0.In order to perform a fair analysis,the CPU times labeled “Uniform”in Table 4 refer to when the less time-consuming,“pre-computed”uniform clouds were used.

We recall that HUA stands for the highest uniform accuracy(achieved with our finest uniform cloud).Let Tγbe the number of seconds spent on attaining the HUA by means of our refinement procedure with a given refinement parameter value γ＞0.

The last column in Table 4 shows the speed-up,Sγ=TU/Tγ.Note that SU=TU/TU=1.When Sγ＞ 1,our refined cloud at level ℓ enables the overall computational cost to be reduced with respect to the uniform discretizationUℓ.

By inspecting the Sγcolumn in Table 4,one can see that high speedups are recorded for Gaussian-based problems P3(GC),P3(GXY5),P3(GXY),P3(GX).

The highest speedups are obtained when dealing with the P3(GX)problem.When γ=10-3,an extremely effective speedup value Sγ＞127.72 is attained.

By inspecting Table 4 one can see that the worst speedups are obtained with the arcTan-based Problem P3(T).Only the row corresponding to γ=5×10-3reports a speedup larger than one.

To explain this not exciting result,we recall that large variations in the solution of problem P3(T)occur over a far broader region than in the other test problems.We recall that in the Gaussian 2D test problems,high variations occur inside a small circle around the Gaussian center,whereas the solution of the “arcTan”2D test problem undergoes high variations on a wider “front”(see Figure 3).Correspondingly,the solutions of the Gaussian 3D test problems undergo non-negligible variations inside a small sphere around the Gaussian center,whereas the solution of problem P3(T)undergoes large variations over a large volume inside[0,1]3,hence many more nodes are needed in order to obtain accurate approximations.

8 Conclusions

An effective,fresh refinement procedure for solving 2D and 3D diffusion problems using DMLPG was introduced and numerically analyzed.

The following points are worth considering.

·Concerning 2D problems:

-our refinement procedure efficiently coalesces nodes around the regions where the solution undergoes large variations.Elsewhere,the discretization is left unchanged.

-A moderate to large improvement in efficiency over uniform clouds is apparent.A far smaller number of nodes may be required to reach a given accuracy by comparison with uniform discretizations.

·As for 3D problems:

-our refinement procedure clusters nodes around the regions of large variation in the solution,as in 2D problems.

-By comparison with uniform discretizations,our refinement procedure allows to significantly reduce the number of nodes,with a consequent saving in CPU running time.

-Let us deal with a test problem in which the Gaussian 3D solution is non-negligible only in the vicinity of the Gaussian “center”.On the other hand,a very high speedup can be achieved by using our refinement procedure instead of uniform clouds.

-Let us assume that the solution of a(test)problem undergoes a large variation on a large portion of the domain.As one can guess,in order to attain a given accuracy,a larger number of nodes is required than in the Gaussian centered case.A smaller gain in efficiency is recorded.

-Linear,tetrahedral FEM requires far more mesh nodes in order to achieve the same accuracy as our DMLPG procedure.As FEM practitioners know,for a given set of discretization nodes,the accuracy of(linear,tetrahedral)FEM strongly depends on the connections in the discretizing mesh.Conversely,the accuracy of our DMLPG-based refinement procedure does not rely on node connections,and it is more robust in response to changes in the discretizations.

Acknowledgement:This work was partially funded by INDAM-GNCS.The third Author received partial support from ADIR funding of Universitá Ca’Foscari Venezia.Some computations were performed on SCSCF,a multiprocessor cluster system owned by Universitá Ca’Foscari Venezia.

Atluri,S.N.(2004):The Meshless method(MLPG)for domain&BIE discretizations.Tech Science,2004,Forsyth GA.

Atluri,S.N.;Zhu,T.(2002): The meshless local Petrov-Galerkin(MLPG)method:A simple&less-costly alternative to the finite element methods.Computer Modeling in Engineering and Sciences,vol.3,no.1,pp.11-51.

Belytschko,T.;Krongauz,Y.;Organ,D.;Fleming,M.;Krysl,P.(1996):Meshless methods:an overview and recent developments.Comp.Methods App.Mech.Eng.,vol.139,pp.3-47.

Chen,L.;Liew,K.(2011):A local Petrov-Galerkin approach with moving Kriging interpolation for solving transient heat conduction problems.Computational Mechanics,vol.47,no.4,pp.455-467.

Fasshauer,G.E.(2007): Meshfree approximation methods with MATLAB,volume 6 of Interdisciplinary Mathematical Sciences.World Scientific Publishing Co.,Singapore.With 1 CD-ROM(Windows,Macintosh and UNIX).

Fries,T.-P.;Matthies,H.-G.(2004): Classification and overview of meshfree methods. Technical Report 2003-3,Technical University Braunschweig,Brunswick,Germany,2004.

Gerace,S.;Erhart,K.;Kassab,A.;Divo,E.(2013): A model-integrated localized collocation meshless method(MIMS). Computer Assisted Methods in Engineering and Science,vol.20,pp.207-225.

Kee,B.B.T.;Liu,G.R.;Zhang,G.Y.;Lu,C.(2008):A residual based error estimator using radial basis functions.Finite Elem.Anal.Des.,vol.44,no.9-10,pp.631-645.

Lancaster,P.;Salkauskas,K.(1981): Surfaces generated by moving least squares methods.Math.Comp.,vol.37,no.155,pp.141-158.

Levin,D.(1998): The approximation power of moving least-squares.Math.Comp.,vol.67,no.234,pp.1517-1531.

Liu,G.R.(2009): Meshfree Methods:Moving Beyond the Finite Element Method.CRC Press,second edition.

Lu,Y.Y.;Belytschko,T.;Gu,L.(1994):A new implementation of the element free Galerkin method.Comp.Methods App.Mech.Eng.,vol.113,pp.397-414.

Mazzia,A.;Pini,G.;Sartoretto,F.(2008): Accurate MLPG solution for 3D potential problems.Computer Modeling in Engineering&Sciences,vol.36,no.1,pp.43-63.

Mazzia,A.;Pini,G.;Sartoretto,F.(2012): Numerical investigation on direct MLPG for 2d and 3d potential problems.Computer Modeling in Engineering&Sciences,vol.88,pp.183-210.

Mazzia,A.;Pini,G.;Sartoretto,F.(2014):Meshless techniques for anisotropic diffusion.APPLIED MATHEMATICS AND COMPUTATION,vol.236C,pp.54-66.

Mazzia,A.;Sartoretto,F.(2010): Meshless solution of potential problems by combining radial basis functions and tensor product ones.Computer Modeling in Engineering&Sciences,vol.68,no.1,pp.95-112.

Mirzaei,D.;Schaback,R.(2013): Direct meshless local Petrov-Galerkin(DMLPG)method:A generalized MLS approximation.Applied Numerical Mathematics,vol.68,no.0,pp.73-82.

Mirzaei,D.;Schaback,R.;Dehghan,M.(2012):On generalized moving least squares and diffuse derivatives.IMA Journal of Numerical Analysis,vol.32,no.3,pp.983-1000.

Strang,G.(2007): Computational Science and Engineering.Wellesley-Cambridge Press,Wellesley,MA.

Yang,S.-W.;Budarapu,P.;Mahapatra,D.;Bordas,S.;Zi,G.;Rabczuk,T.(2015): A meshless adaptive multiscale method for fracture. Computational Materials Science,vol.96,no.PB,pp.382-395.cited By 1.

1Dipartimento ICEA,Università di Padova,Via Trieste 63,35121 Padova,Italy{annamaria.mazzia,giorgio.pini}@unipd.it

2DAIS,Università Ca’Foscari Venezia,Via Torino 155,10173 Venezia,Italy flavio.sartoretto@unive.it