A generic, electro-thermo-mechanically coupled finite element program is developed for three-dimensional simulation of resistance welding. The developed computer program has reached a level of a complete standalone software that can be utilized as a tool in the analysis of resistance welding processes. Contact between deformable objects is modeled by the penalty method to handle multiple objects that appear in joining processes. Two algorithms are implemented for the identification of contact pairs.Simulations are presented for two metal forming processes involving contact and a number of resistance welding processes, which cover a wide range of spot welding and projection welding applications. Three-dimensional simulation of spot welding enables the analysis of critical effects like electrode misalignment and shunt effects between consecutive spots. A single-sided spot welding case involving three-dimensional contact is also presented. This case was suggested by and discussed with a German steel manufacturer. When it comes to projection welding, a natural need for three-dimensional analysis arises in many cases because of the involved geometries. Cross-wire welding and welding of square nuts to sheets by projection welding are presented by means of experiments and simulations. These two cases are used to explore the capabilities of the developed simulation software by comparing experiments and simulations. A number of other projection welding cases are presented for further application of the software. These include joining of parallel sheets by circular projections and joining of perpendicular sheets by longitudinal projections. In the former case, the effects of unequal projection heights are analyzed, and in the latter case, the simulations are compared to the corresponding experiment by a Japanese company that proposed the case. Another industrial case, by a German company, is joining of micro components. The joining is based on mechanical locking, and the deformation is accommodated by resistance heating, which at the same time is used to melt a polymer coating locally for creating electrical contact necessary for the end-product. All the above cases are modeled by meshing techniques included in the computer program. Structured, isoparametric meshing is utilized for setting up initial meshes of individual objects. Unstructured, all-hexahedral meshing is utilized for creating initial meshes of objects defined by CAD surfaces and is applied for remeshing of selected objects in order to carry on certain simulations. The all-hexahedral meshing procedures are enhanced by adaptive bounding boxes, facilities for handling multi-object simulations and overall improved by applying topology based criteria in the creation of hexahedral meshes. Simulation time is significantly reduced by a developed parallel skyline solver. The new solver is developed for shared memory and can be implemented in existing finite element codes by changing the call to the solver, as long as the system matrix is prepared in skyline format. Finally, the above models and procedures are operated by a developed graphical interface including its own pre and post processing facilities. This combines the above into a new complete, standalone software: SORPAS 3D.