Coronary blood flow is different from the flow in the other parts of the cardiovascular system because coronary flow is strongly influenced by the contraction and relaxation of the ventricles and atria. The ventricular and atrial pressure can be used to approximate the intramyocardial pressure and model the compressive force acting on the coronary vessels caused by the contraction and relaxation of the ventricles and atria throughout the cardiac cycle. The ventricular and atrial pressure should be determined by the interactions between the heart and the vascular system. Most of the previous studies on three-dimensional computations of coronary flow ignored the interactions between the heart and the vascular system and mainly solved for flow only in the coronary arteries. Prior studies that considered the interactions between the heart and the vascular impedance, and computed coronary flow by considering the ventricular pressure as the intramyocardial pressure, were applied to idealized, rigid three-dimensional models with low mesh resolutions. Additionally, the interaction was considered using an explicit coupling approach, which requires either sub-iterations within the same time step or an explicit time integration scheme that bounds a time step size for stability. To compute physiologically realistic flow and pressure, a robust and stable flow solver that models aortic pressure accurately and couples coronary impedance and intramyocardial pressure at the coronary outlet boundaries is needed. In this study, I developed two sets of new boundary conditions with constraints on the shape of the velocity profiles to satisfy these needs. The first boundary condition couples a system composed of two lumped parameter heart and venous models and a lumped parameter pulmonary model to a three-dimensional finite element model of the aorta with appropriate outlet boundary conditions. One of the two lumped parameter heart models, which models the left side of the heart, is coupled to the inlet of a three-dimensional finite element model of the aorta to model the interactions between the heart and the systemic circulation. The second boundary condition couples lumped parameter coronary vascular models to coronary outlets of the three-dimensional finite element model. The ventricular and atrial pressure computed with the two lumped parameter heart models can be used to represent the intramyocardial pressure of the coronary vascular beds. Additionally, an augmented Lagrangian method was implemented to enforce constraints on the shape of the velocity profiles on the Neumann boundaries with complex flow structures to obtain robust and stable solutions. Both of these boundary conditions together with the constraints on the velocity profiles enable us to compute physiologic coronary flow and pressure, aortic flow and pressure, and ventricular pressure and volume. In this thesis, I presented an additional application of this method that models autoregulatory mechanisms of the cardiovascular system. This was accomplished by designing control loop systems that couple the lumped parameter heart models, arterial system, and coronary vascular models. This method can be used to study the interactions between the heart, coronary vascular beds, and arterial system and investigate improved treatments for cardiovascular disease.