Disclosed are techniques for representing and modeling one or more systems in which each system corresponds to an application mode. This may be done for one or more geometries using local and/or non-local couplings. For each application mode, physical quantities are modeled and may be defined using a graphical user interface. Physical properties may be used to model the physical quantities of each system. The physical properties may be defined in terms of numerical values or constants, and mathematical expressions that may include numerical values, space coordinates, time coordinates, and actual physical quantities. Physical quantities and any associated variables may apply to some or all of a geometric domain, and may also be disabled in other parts of a geometrical domain. Partial differential equations describe the physical quantities. One or more application modes may be combined using an automated technique into a combined system of partial differential equations as a multiphysics model. A portion of the physical quantities and variables associated with the combined system may be selectively solved for. Also described are methods for computing the stiffness matrix, residual vector, constraint matrix, and constraint residual vector for the finite element discretization of a system of partial differential equations in weak form that includes local and non-local variables coupling multiple geometries.