Modern manufacturing engineering is based on a ``design-through-analysis'' workflow. According to this paradigm, a prototype is first designed with Computer-aided-design (CAD) software and then finalized by simulating its physical behavior, which usually involves the simulation of Partial Differential Equations (PDEs) on the designed product. The simulation of PDEs is often performed via finite element discretization techniques. A severe bottleneck in the entire process is undoubtedly the interaction between the design and analysis phases. The prototyped geometries must undergo the time-consuming and human-involved meshing and feature removal processes to become ``analysis-suitable''. This dissertation aims to develop and study numerical solvers for PDEs to improve the integration between numerical simulation and geometric modeling. The thesis is made of two parts. In the first one, we focus our attention on the analysis of isogeometric methods which are robust in geometries constructed using Boolean operations. We consider geometries obtained via trimming (or set difference) and union of multiple overlapping spline patches. As differential model problems, we consider both elliptic (the Poisson problem, in particular) and saddle point problems (the Stokes problem, in particular). As it is standard, the Nitsche method is used for the weak imposition of the essential boundary conditions and to weakly enforce the transmission conditions at the interfaces between the patches. After proving through well-constructed examples that the Nitsche method is not uniformly stable, we design a minimal stabilization technique based on a stabilized computation of normal fluxes (and on a simple modification of the pressure space in the case of the Stokes problem). The main core of this thesis is devoted to the derivation and rigorous mathematical analysis of a stabilization procedure to recover the well-posedness of the discretized problems independently of the geometric configuration in which the domain has been constructed. In the second part of the thesis, we consider a different approach. Instead of considering the underlying spline parameterization of the geometrical object, we immerse it in a much simpler and readily meshed domain. From the mathematical point of view, this approach is closely related to the isogeometric discretizations in trimmed domains treated in the first part. In this case, we consider the Raviart-Thomas finite element discretization of the Darcy flow. First, we analyze a Nitsche and a penalty method for the weak imposition of the essential boundary conditions on a boundary fitted mesh, a problem that was not studied before, not needed for our final goal, but still interesting by itself. Then, we consider the case of a general domain immersed in an underlying mesh unfitted with the boundary. We focus on the Nitsche method presented for the boundary fitted case and study its extension to the unfitted setting. We show that the so-called ghost penalty stabilization provides an effective solution to recover the well-posedness of the formulation and the well-conditioning of the resulting linear system.