Usually the FDM is a second order method on a rectangular grid. In the FIDISOL (Finite Difference Solver) program package for the solution of nonlinear systems of elliptic and parabolic PDEs (black-box), see [1], Chapter 17, we have generalized the FDM to arbitrary consistency order. The access to the discretization error is by the difference of different consistency orders. We derive an error equation that makes transparent the influence of all types of discretization and linearization errors and is the key to selfadaptation of the whole solution process, e.g. to select the optimal consistency order, space and time mesh size, stopping criterion for Newton iteration and iterative linear solver. In order to increase the geometrical flexibility we developed the CADSOL (Cartesian Arbitrary Domain Solver) package that solves the PDEs on an arbitrary domain with body-oriented grid (rectangular in index space). Here we learned how to generate difference formulas of arbitrary order on an arbitrary set of nodes by influence functions and we introduced dividing lines (DLs) that separate subdomains with different PDEs and couple the solution globally by coupling conditions (CCs). However, the body-oriented grid is a severe restriction. In the FDEM (Finite Difference Element Method) program package, see [2], we generalize our solution method to an arbitrary 2-D/3-D domain with triangular/tetrahedral grid. The consistency order is arbitrary, the FEM grid gives only the space structure. A sophisticated algorithm has been developed to select nodes in rings (2-D) or balls (3-D) for the generation of the difference formulas. The same error equation allows mesh refinement and order selection. By DLs and CCs we treat a domain composed of subdomains with different PDEs.

Finally in the S-FDEM (Sliding FDEM) package (under development) we generalize our solution method to a domain whose subdomains may slide relative to the other subdomains, separated by S-DLs (sliding DLs). In the different domains are quite different meshes, so one has non-matching moving grids, with individual mesh refinement in the subdomains and global error estimate. These requests result from industry cooperations, e.g. if there is for a high pressure diesel injection pump the large housing, the thin lubrication gap and the rotating shaft, the whole can be solved uniformly with global error estimate (note that we have a black-box solver). The (S-)FDEM program package is parallelized with optimal data structures for distributed memory parallel computers, see [3]. The (iterative) solution of the resulting large and sparse linear systems is made by the LINSOL program package, see [4], that has polyalgorithms with automatic method switching between generalized conjugate gradient methods of different properties. This is even a selfadaptation of the solution method.