Accuracy of numerical integration in material point method-based geotechnical analysis and simulation by optimizing integration weights

Abstract

In one embodiment, a technique for numerical integration in material point method (MPM)-based geotechnical analysis and simulation is provided. Input terms for an element of a background mesh are received. The input terms including material points in the element that describe a continuum of soil, rock and/or groundwater. A set of constraints is created that defines an optimization problem. The set of constraints provide that numerical integration of the material points weighted by unknown integration weights equal exact integration for finite element shape functions. The optimization problem defined by the constraints is solved by an optimization algorithm to minimize numerical integration error for polynomials up to a given order to produce a set of integration weights. The set of integration weights is scaled to conserve the mass of the material points to produce optimized integration weights. The optimized integration weights are used in numerical integration performed in MPM-based geotechnical analysis and simulation.

Claims

1. A method for numerical integration in material point method (MPM)-based geotechnical analysis and simulation, comprising: receiving, by geotechnical engineering software executing on one or more computing devices, input terms for an element of a background mesh, the input terms including material points in the element that describe a continuum of soil, rock and/or groundwater; conducting, by the geotechnical engineering software, a MPM-based geotechnical analysis and simulation of the background mesh that includes the element at least in part by: creating a set of constraints that provide that, for a plurality of finite element shape functions up to a given order, numerical integration of the plurality of finite element shape functions at the material points weighted by unknown integration weights equal exact integration for the plurality of finite element shape functions over an area of the element, wherein the given order is an order greater than one, performing an optimization, by an optimization algorithm of the geotechnical engineering software, based on the set of constraints to determine the unknown integration weights, the optimization algorithm to produce a set of integration weights by minimizing numerical integration error for polynomials up to the given order, scaling the set of integration weights to conserve the mass of the material points to produce optimized integration weights, and performing a given order numerical integration at the material points weighted by the optimized integration weights; and displaying in a user interface, or storing to a computing device-readable medium, results of the MPM-based geotechnical analysis and simulation.

2. The method of claim 1, wherein the background mesh includes a plurality of elements, and at least the receiving, creating, performing the optimization and scaling are performed for each of the plurality of elements.

3. The method of claim 1, wherein the input terms further include a signed distance field whose value depends on location outside or inside a surface of a point cloud defined by the material points, and the set of constraints provide that an integration domain of the numerical integration is a region inside or on the surface of the point cloud.

4. The method of claim 1, wherein the input terms further include an original total mass of the material points in the element, and the scaling comprises applying a difference between the original total mass and current total mass to the set of integration weights.

5. The method of claim 1, wherein the set of constraints include constraints: ${.Math.}_{i=1}^n{w}_i{}_0\left(p_i\right)={}_{}{}_{0^{-}})$ $.Math.$ ${.Math.}_{i=1}^n{w}_i{}_m\left(p_i\right)={}_{}{}_{m^{-}})$ where n is a number of material points in the element, w.sub.1, . . . , w.sub.n are integration weights, for material points, .sub.0 . . . .sub.m are finite element shape functions, p.sub.1, . . . , p.sub.n are material points and is an integration domain equal to the extent of the element.

6. The method of claim 1, wherein the optimization algorithm is a least squares algorithm that minimizes quadratic error in the Euclidean norm.

7. The method of claim 1, wherein the optimization algorithm constrains each integration weight of the set of integration weights to be non-negative.

8. The method of claim 1, wherein at least some of the integration weights of the set of optimized integration weights differ from each other.

9. The method of claim 1, wherein the element is only partially covered by material points.

10. The method of claim 1, wherein the performing the given order numerical integration applies the optimized integration weights directly to the material points.

11. A computing device comprising: a display screen; a processor; and a memory coupled to the processor and configured to store geotechnical engineering software, the geotechnical engineering software including a mesher software module configured to generate a background mesh composed of elements formed from nodes, each element including material points that describe a continuum of soil, rock and/or groundwater, a solver software module configured to use material point method (MPM) to analyze and simulate behavior of the soil, rock and/or groundwater at least in part by performing a given order numerical integration at the material points on each element, the numerical integration weighted by a set of optimized integration weights, and a numerical integration software sub-module configured to determine the set of optimized integration weights used to perform the numerical integration by creating a set of constraints that provide that, for a plurality of finite element shape functions up to the given order, numerical integration of the plurality of finite element shape functions at the material points weighted by unknown integration weights equal exact integration for the plurality of finite element shape functions over an area of a respective element, wherein the given order is an order greater than one, performing an optimization based on the set of the constraints with an optimization algorithm to determine the unknown integration weights, the optimization to produce a set of integration weights by minimizing numerical integration error of polynomials up to the given order, and scaling the set of integration weights to conserve the mass of the material points to produce the set of optimized integration weights.

12. The computing device of claim 11, wherein the set of constraints provide that an integration domain of the numerical integration is a region inside or on a surface of a point cloud defined by the material points.

13. The computing device of claim 11, wherein the scaling applies a difference between an original total mass and current total mass to the set of integration weights.

14. A non-transitory computing device readable medium having instructions stored thereon, the instructions when executed by one or more computing devices operable to: receive input terms for an element of a background mesh, the input terms including material points in the element that describe a continuum of soil, rock and/or groundwater; conduct a material point method (MPM)-based geotechnical analysis and simulation of the background mesh that includes the element at least in part by: creating a set of constraints that provide that, for a plurality of finite element shape functions up to a given order, numerical integration of the plurality of finite element shape functions at the material points weighted by unknown integration weights equal exact integration for the plurality of finite element shape functions over an area of the element, wherein the given order is an order greater than one, performing an optimization based on the constraints with an optimization algorithm to determine the unknown integration weights, the optimization algorithm to produce a set of integration weights by minimizing numerical integration error for polynomials up to the given order, scaling the set of integration weights to conserve the mass of the material points to produce optimized integration weights, and performing a given order numerical integration at the material points weighted by the optimized integration weights; and display in a user interface, or store to a computing device-readable medium, results of the MPM-based geotechnical analysis and simulation.

15. The non-transitory electronic-device readable medium of claim 14, wherein the background mesh includes a plurality of elements, and at least the instructions that when executed are operable to receive and conduct are operable to perform operations for each of the plurality of elements.

16. The non-transitory electronic-device readable medium of claim 14, wherein the input terms further include a signed distance field whose value depends on location outside or inside a surface of a point cloud defined by the material points, and the set of constraints provide that an integration domain of the numerical integration is a region inside or on the surface of the point cloud.

17. The non-transitory electronic-device readable medium of claim 14, wherein the input terms further include an original total mass of the material points in the element, and the instructions that when executed are operable to scale include instructions that when executed are operable to apply a difference between the original total mass and current total mass to the set of integration weights.

18. The non-transitory electronic-device readable medium of claim 14, wherein the set of constraints include constraints: ${.Math.}_{i=1}^n{w}_i{}_0\left(p_i\right)={}_{}{}_{0^{-}})$ $.Math.$ ${.Math.}_{i=1}^n{w}_i{}_m\left(p_i\right)={}_{}{}_{m^{-}})$ where n is a number of material points in the element, w.sub.1, . . . , w.sub.n are integration weights, for material points, .sub.0 . . . .sub.m are finite element shape functions, p.sub.1, . . . , p.sub.n are material points and is an integration domain equal to the extent of the element.

19. The non-transitory electronic-device readable medium of claim 14, wherein the optimization algorithm is a least squares algorithm that minimizes quadratic error in the Euclidean norm.

20. The non-transitory electronic-device readable medium of claim 14, wherein the optimization algorithm constrains each integration weight of the set of integration weights to be non-negative.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The description below refers to the accompanying drawings of example embodiments, of which:

(2) FIG. 1 is a high-level block diagram of an example software architecture for geotechnical engineering software that may perform MPM-based geotechnical analysis and simulation;

(3) FIG. 2 is a diagram of an example computational mesh (background mesh) and material points; and

(4) FIG. 3 is a flow diagram of an example sequence of steps that may be performed to improve accuracy of numerical integration by optimizing integration weights.

DETAILED DESCRIPTION

(5) FIG. 1 is a high-level block diagram of an example software architecture 100 for geotechnical engineering software that may perform MPM-based geotechnical analysis and simulation. The geotechnical engineering software may be a stand-alone software application or a component of a larger software application. In one example implementation, the geotechnical engineering software is the PLAXIS3D geotechnical engineering software available from Bentley Systems of Exton, PA. However, it should be understood that the geotechnical engineering software may take a variety of other forms.

(6) The geotechnical engineering software may be divided into local software 110 that executes on one or more computing devices local to an end-user (collectively local devices) and, in some cases, cloud-based software 120 that is executed on one or more computing devices remote from the end-user (collectively cloud computing devices) accessible via a network (e.g., the Internet). Each computing device may include processors, memory/storage, a display screen, and other hardware (not shown) for executing software, storing data and/or displaying information. The local software 110 may include a number of software modules operating on a local device and the cloud-based software 120, if present, may include, additional software modules operating on cloud computing devices. Tasks may be divided in a variety of different manners among the software modules. For example, in one implementation, software modules of the local software 110 may be responsible for performing non-processing intensive operations and software modules of the cloud-based software 120 may perform more processing intensive operations. However, many other arrangements may be employed.

(7) The software modules may include a user interface module 130, a mesher module 140 and a solver module 150. The user interface module 130 may be responsible for providing a user interface (e.g., a graphical user interface and/or a command line interface,) for receiving user input and displaying output. The mesher module 140 may be responsible for generating a computational mesh (background mesh) composed of elements formed from nodes that model a continuum of soil, rock and/or groundwater, for example based on information (e.g., a model, geometric description, parameters/conditions, etc.) selected or provided in the user interface. The solver module 150 may be responsible for using MPM to analyze and simulate behavior (e.g., deformation, stability, interactions, and the like) of the soil, rock and/or groundwater. To that end, the solver module 150 may associate material points with information (e.g., mass, volume, stress, state variables, etc.) and move the material points through elements of the background mesh over time steps of analysis and simulation. The solver module may include a numerical integration sub-module 155 that is used in these operations. Results of the analysis and simulation may be displayed by the user interface module 130, stored in a file, database, etc., provided to other software, sand the like.

(8) FIG. 2 is a diagram 200 of an example computational mesh (background mesh) 210 and material points 220 that may assist in understanding the operation of the mesher module 140 and solver module 150. The background mesh 210 and material points 220 may represent a continuum of soil, rock and/or groundwater. The background mesh 210 is composed of elements (e.g., triangles) formed from nodes. The material points 220 are distributed within the elements. Some elements 210 may be fully covered by material points, while other elements may be only partially covered.

(9) When performing MPM-based geotechnical analysis and simulation using a background mesh and material points such as shown in the example in FIG. 2, the solver module 150, working together with a numerical integration sub-module 155, may perform numerical integration over each element at each time step of the analysis or simulation to obtain entries of a system matrix and right-hand side vector. An example integration may be given as:
.sub.KN.sub.i.sub.pw.sub.pN.sub.i(x.sub.p)
where K is an element of the background mesh, N.sub.i is a finite element shape function, x.sub.p is the position of a material point p, and w.sub.p is the integration weight of the material point p.

(10) In one or more embodiments, the accuracy of such a numerical integration used is improved by using optimized integration weights. The original material points may be used for integration, rather than inter-extrapolating information from material points to nodes of the elements of the background mesh. Differing weights for each material point may be used, rather than using all equal weights.

(11) FIG. 3 is a flow diagram of an example sequence of steps 300 that may be performed by the geotechnical engineering software to improve accuracy of numerical integration by optimizing integration weights. The sequence of step 300 are described in relation to a single element, however it should be understood that such steps may be executed in a loop (or in parallel) over all non-empty elements of a background mesh. At step 310, the numerical integration sub-module 155 receives input terms (e.g., from the solver module 150) for the element. The input terms may include positions of the material points x.sub.1, x.sub.2, . . . , x.sub.n in the element, a signed distance field L, and an original total mass mx of the material points in the element. The signed distance field L is a continuous function that indicates the distance of a given point x from the boundary of a point cloud defined by the material points, with the sign determined by whether it is outside or inside the point cloud (e.g., positive outside and negative inside). The zero iso-surface of the signed distance field describes the surface of the point cloud. In one implementation, the signed distance field may be a level-set field. However, it should be understood that it may also take other forms.

(12) At step 320, the numerical integration sub-module 155 creates a set of constraints that define an optimization problem. Since the number of constraints usually does not coincide with the number of integration weights to be determined, the resulting system generally cannot be solved directly, hence it is considered an optimization problem. In general, the constraints may provide that numerical integration of the material points weighted by unknown integration weights equal exact integration for finite element shape functions.

(13) To further understand the operation of step 320, the numerical integration may be examined in more detail. As discussed above, in numerical integration a sum replaces an integral such that the integrand is sampled at a number of discrete integration points. As such, if a function is integrated over an integration domain 22, numerical integration may be represented as:
f.sub.pw.sub.pf(x.sub.p)
where {x.sub.p} are the positions of integration points (here, the positions of material points) and {w.sub.p} are the integration weights. Basically, the function is evaluated at the integration points (material points) and the obtained function values are weight-accumulated with respect to the integration weights. The integration weights may depend on the volume of the integration domain 2 as well as the chosen numerical integration method. Typically, integration weights are different for each integration point in an accurate higher-order numerical integration method (e.g., Gauss quadrature).

(14) When numerical integration is used to obtain the entries of the system matrix and right-hand side vector, the integration domain 2 may be chosen to be the region (e.g., area, volume, etc.) of element K that is covered by the material point cloud (i.e. inside or on the surface of the material point cloud). In other words, the part of K in which the signed distance function L is negative, may be represented as:
={xK:L(x)0}.
A set of constraints may be constructed that ensures that polynomials up to order k can be integrated accurately in the element.

(15) To construct the set of constraints, function may be chosen from finite element shape functions N.sub.0, . . . , N.sub.m up to order k, such that there are m+1 possible finite element shape functions. The finite element shape functions may be generated as standard shape functions or hierarchical shape functions. In a simple example, shape functions in various orders may be generated by products of appropriate polynomials of the coordinate system of the element. In other examples, shape functions may be generated to be dependent on nodal values placed on the element boundary. It should be understood that a wide variety of known techniques may be used to generate the finite element shape functions, and that the present disclosure is not limited to use with one family or category of shape function.

(16) With m+1 possible finite element shape functions, in one implementation, a set of constraints may be defined as:

(17) ${.Math.}_{i=1}^n{w}_i{}_0\left(p_i\right)={}_{}{}_{0^{-}})$$.Math.$${.Math.}_{i=1}^n{w}_i{}_m\left(p_i\right)={}_{}{}_{m^{-}}).$
On the right-hand side, we have an exact integration of finite element shape functions .sub.0 . . . .sub.m over a domain equal to the extent (e.g., area) of the element. The exact integration may be obtained, for example, by a very dense subdivided Gau integration, or another integration technique. On the left-hand side we have a numerical integration of the shape functions .sub.0, . . . , .sub.m using the provided set of material points p.sub.1, . . . , p.sub.n and yet unknown integration weights w.sub.1, . . . , w.sub.n for each of n material points

(18) Such a set of constraints may be written in a matrix-vector notation as Aw=b, where:

(19) $A=\left(\begin{array}{ccc}{}_0\left(p_1\right)& .Math.& {}_0\left(p_n\right)\\ .Math.& & .Math.\\ {}_m\left(p_0\right)& .Math.& {}_m\left(p_n\right)\end{array}\right)$
is a m+1 by n+1 matrix, and:

(20) $b=\left(\begin{array}{cc}{}_{}& {}_0\\ {}_{}& {}_0\\ {}_{}& {}_m\end{array}\right)$
is a right-hand side vector, and w is a sought vector of integration weights.

(21) While the specific definitions above may provide constraints that ensure numerical integration of the material points weighted by unknown integration weights equal exact integration for finite element shape functions, it should be understood that in other implementations other specific definitions may be used, and that the techniques are not limited to these specific definitions.

(22) At step 330, the numerical integration sub-module 155 uses an optimization algorithm to solve the optimization problem created in step 320. The optimization algorithm may minimize numerical integration error for polynomials up to a given order to produce a set of integration weights. A variety of different optimization algorithms may be employed. In one implementation, the optimization algorithm is a least-squares algorithm that attempts to minimize the quadratic error in the Euclidean norm (the square root of the sum of the squares of the coordinates) given as:

(23) $\underset{w}{\min }{.Math.\mathrm{Aw}-b.Math.}_2^2$
where, A is the m+1 by n+1 matrix, w is the vector of integration weights, and b is the right-hand side vector.

(24) In general, it may be desirable that integration weights w.sub.1, . . . , w.sub.n be non-negative (i.e. w0) to ensure numerical stability of the integration method. Therefore, the above optimization problem may be modified to also include this constraint, and thus be given as:

(25) $\underset{w0}{\min }{.Math.\mathrm{Aw}-b.Math.}_2^2$
The optimization problem may be solved using a least-squares algorithm (e.g., using a conventional non-negative least-squares solver) to obtain a set of integration weights w.sub.1, . . . , w.sub.n that are non-negative and integrate all polynomials of up to degree k on integration domain as accurately as possible with the given material points x.sub.1, x.sub.2, . . . , x.sub.n subject to integration.

(26) While a least-squares optimization above may be used to solve the optimization problem created in step 320, it should be understood that other optimization algorithms may be used to solve the problem. For example, the optimization algorithm may be a min-max (Chebyshev) algorithm, a minimum sum algorithm, or another type of optimization algorithm. Such alternative algorithms may use alternative formulations differing from the specific examples shown above.

(27) At step 340, the numerical integration sub-module 155 scales the integration weights to conserve the mass of the material points to produce optimized integration weights. Scaling of the integration weights may ensure that mass of the material points x.sub.1, x.sub.2, . . . , x.sub.n in element K are preserved. In one implementation, this may be achieved by scaling the optimized integration weights {w.sub.p} determined in step 330 with a correction factor. Current total mass m of the material points in the element may be given as:
m=.sub.pw.sub.p.
A correction factor may be given as

(28) $\frac{m_K}{m}$
where m.sub.K is the original total mass of the material points in the element. The correction factor may be applied to the integration weights (i.e. multiplied against them) to yield optimized integration weights that conserve mass.

(29) While a specific correction factor formulation is provided above, it should be understood that other formulations may be used to produce optimized integration weights that conserve mass.

(30) At step 350, the optimized integration weights determined in step 340 are used by the solver 150 of the geotechnical engineering software in a numerical integration performed as part of MPM-based geotechnical analysis and simulation. The results of such MPM-based geotechnical analysis and simulation may be displayed by the geotechnical engineering software in a user interface, stored to an computing device-readable medium, or otherwise used.

(31) It should be understood that various adaptations and modifications may be readily made to what is described above, to suit various implementations and environments. While it is discussed above that certain steps may be implemented in a specific sequence, it should be understood that other sequences may be used, and that such sequences may include additional, intervening steps, steps implemented in parallel, and the like. While it is discussed above that aspects of the techniques may be implemented by specific software processes executing on specific hardware, it should be understood that some or all of the techniques may also be implemented by different software on different hardware. In addition to general-purpose computing devices, the hardware may include specially configured logic circuits and/or other types of hardware components. Above all, it should be understood that the above descriptions are meant to be taken only by way of example.