Physics-based simulation, coupled with the Finite Element Method (FEM), has emerged asa powerful tool in understanding and analyzing the complex behavior of elastic materials. Implicit Discretization is essential for efficiently and accurately simulating elastic solids and cloth. In this thesis, we first explore ways of creating volumetric mesh for embedding surface mesh. The embedded surface mesh has many small self intersections. We devise an efficient and robust way of generating a hexahedron mesh to embed the triangle mesh, so that the self-intersecting regions are correctly duplicated. We then simulate the hexahedron mesh using Finite Element Method and interpolate to the embedded triangle mesh. The second part explores different ways of improving the core simulation solver of Finite Element Method. We improve on the Position Based Dynamics (PBD)/Extended Position Based Dynamics (XPBD) framework. PBD/XPBD are known for their robustness under a very small computational budget. However, they have several limitations. PBD/XPBD does not converge when the computational budget is sufficient. PBD/XPBD also supports limited hyperelastic constitutive models. We devise PXPBD (Primary Extended Position Based Dynamics) to address these issues. The PXPBD methods improves convergence rate of PBD/XPBD and support arbitrary hyperelastic models. PBD/XPBD framework does not support quasistatic simulation either. We note quasistatic simulation is important for generating training data for machine learning based simulation approaches such as QNN (Quasistatic Neural Network). We also notice that the constraint-based Gauss-Seidel approach in PBD/XPBD causes loss of information on the node solve. So we design Position-Based Nonlinear Gauss Seidel (PBNG), where the hyperelastic energy from FEM is optimized per node, instead of per constraint. Doing so not only boosts convergence rate, but also enables high quality quasistatic simulation, which runs much faster than the existing methods such as Newton’s method. The last part of my thesis uses a machine learning-based approach to simulate human musculature effectively. We use a neural network to model the deformation of human soft tissues such as muscle, tendon, fat and skin with high fidelity. By deploying a biomechanicsbased approach, we estimate the activation of single muscle during deformation. The activation parameters are then incorporated into an active neural network (ANN). which adds per-vertex deformation on top of Linear Blend Skinning (LBS). Our neural network achieves more than 1000X speed up as compared to traditional FEM approaches, but it has the same level of visual fidelity.