diff --git a/NAMESPACE b/NAMESPACE
index f3871ec..fc48101 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -25,6 +25,7 @@ export(layout_with_centrality_group)
 export(layout_with_constrained_stress)
 export(layout_with_constrained_stress3D)
 export(layout_with_eigen)
+export(layout_with_fixed_coords)
 export(layout_with_focus)
 export(layout_with_focus_group)
 export(layout_with_pmds)
diff --git a/NEWS.md b/NEWS.md
index e895f58..ce41cc2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,7 @@
   work for disconnected graphs
 * internal code refactoring
 * added `layout_as_metromap()`
+* added `layout_with_fixed_coords()`
 
 # graphlayouts 1.0.2
 
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 7f4a637..50ca6a4 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -9,6 +9,10 @@ constrained_stress_major <- function(y, dim, W, D, iter, tol) {
     .Call(`_graphlayouts_constrained_stress_major`, y, dim, W, D, iter, tol)
 }
 
+fixed_stress_major <- function(y, fixedCoords, W, D, iter, tol) {
+    .Call(`_graphlayouts_fixed_stress_major`, y, fixedCoords, W, D, iter, tol)
+}
+
 constrained_stress3D <- function(x, W, D) {
     .Call(`_graphlayouts_constrained_stress3D`, x, W, D)
 }
diff --git a/R/layout_stress.R b/R/layout_stress.R
index 3bbb72a..1c3f40c 100644
--- a/R/layout_stress.R
+++ b/R/layout_stress.R
@@ -438,6 +438,77 @@ layout_with_constrained_stress3D <- function(g, coord, fixdim = "x", weights = N
         ))
     }
 }
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#' Layout with fixed coordinates
+#'
+#' @name layout_fixed_coords
+#' @description force-directed graph layout based on stress majorization with
+#' fixed coordinates for some nodes
+#' @param g igraph object
+#' @param coords numeric n x 2 matrix, where n is the number of nodes. values
+#' are either NA or fixed coordinates. coordinates are only calculated for the
+#' NA values.
+#' @param weights possibly a numeric vector with edge weights. If this is NULL and the graph has a weight edge attribute, then the attribute is used. If this is NA then no weights are used (even if the graph has a weight attribute). By default, weights are ignored. See details for more.
+#' @param iter number of iterations during stress optimization
+#' @param tol stopping criterion for stress optimization
+#' @param mds should an MDS layout be used as initial layout (default: TRUE)
+#' @param bbox constrain dimension of output. Only relevant to determine the placement of disconnected graphs
+#' @details Be careful when using weights. In most cases, the inverse of the edge weights should be used to ensure that the endpoints of an edges with higher weights are closer together (weights=1/E(g)$weight).
+#'
+#' The layout_igraph_* function should not be used directly. It is only used as an argument for plotting with 'igraph'.
+#' 'ggraph' natively supports the layout.
+#' @seealso [layout_constrained_stress]
+#' @return matrix of xy coordinates
+#' @examples
+#' library(igraph)
+#' set.seed(12)
+#' g <- sample_bipartite(10, 5, "gnp", 0.5)
+#' fxy <- cbind(c(rep(0, 10), rep(1, 5)), NA)
+#' xy <- layout_with_fixed_coords(g, fxy)
+#' @export
+layout_with_fixed_coords <- function(g, coords, weights = NA,
+                                     iter = 500, tol = 0.0001, mds = TRUE, bbox = 30) {
+    ensure_igraph(g)
+
+    oldseed <- get_seed()
+    set.seed(42) # stress is deterministic and produces same result up to translation. This keeps the layout fixed
+    on.exit(restore_seed(oldseed))
+
+    if (missing(coords)) {
+        stop('"coords" is missing with no default.')
+    }
+    if (nrow(coords) != igraph::vcount(g) && ncol(coords) != 2) {
+        stop("coords has the wrong dimensions")
+    }
+    if (all(!is.na(coords))) {
+        warning("all coordinates fixed")
+        return(coords)
+    }
+    comps <- igraph::components(g, "weak")
+    if (comps$no == 1) {
+        if (igraph::vcount(g) == 1) {
+            return(matrix(c(0, 0), 1, 2))
+        } else {
+            D <- igraph::distances(g, weights = weights)
+            W <- 1 / D^2
+            diag(W) <- 0
+            n <- igraph::vcount(g)
+            xinit <- .init_layout(g, D, mds, n, dim = 2)
+            not_na_idx <- which(!is.na(coords), arr.ind = TRUE)
+            xinit[not_na_idx] <- coords[not_na_idx]
+            return(fixed_stress_major(xinit, coords, W, D, iter, tol))
+        }
+    } else {
+        # return(.component_layouter(
+        #     g = g, weights = weights, comps = comps, dim = 2, mds = mds,
+        #     bbox = bbox, iter = iter, tol = tol, FUN = constrained_stress_major, fixdim = fixdim, coord = coord
+        # ))
+        stop("only connected graphs are supported")
+    }
+}
+
 #' radial focus group layout
 #'
 #' @description arrange nodes in concentric circles around a focal node according to their distance from the focus and keep predefined groups in the same angle range.
diff --git a/man/layout_fixed_coords.Rd b/man/layout_fixed_coords.Rd
new file mode 100644
index 0000000..5734e5a
--- /dev/null
+++ b/man/layout_fixed_coords.Rd
@@ -0,0 +1,57 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/layout_stress.R
+\name{layout_fixed_coords}
+\alias{layout_fixed_coords}
+\alias{layout_with_fixed_coords}
+\title{Layout with fixed coordinates}
+\usage{
+layout_with_fixed_coords(
+  g,
+  coords,
+  weights = NA,
+  iter = 500,
+  tol = 1e-04,
+  mds = TRUE,
+  bbox = 30
+)
+}
+\arguments{
+\item{g}{igraph object}
+
+\item{coords}{numeric n x 2 matrix, where n is the number of nodes. values
+are either NA or fixed coordinates. coordinates are only calculated for the
+NA values.}
+
+\item{weights}{possibly a numeric vector with edge weights. If this is NULL and the graph has a weight edge attribute, then the attribute is used. If this is NA then no weights are used (even if the graph has a weight attribute). By default, weights are ignored. See details for more.}
+
+\item{iter}{number of iterations during stress optimization}
+
+\item{tol}{stopping criterion for stress optimization}
+
+\item{mds}{should an MDS layout be used as initial layout (default: TRUE)}
+
+\item{bbox}{constrain dimension of output. Only relevant to determine the placement of disconnected graphs}
+}
+\value{
+matrix of xy coordinates
+}
+\description{
+force-directed graph layout based on stress majorization with
+fixed coordinates for some nodes
+}
+\details{
+Be careful when using weights. In most cases, the inverse of the edge weights should be used to ensure that the endpoints of an edges with higher weights are closer together (weights=1/E(g)$weight).
+
+The layout_igraph_* function should not be used directly. It is only used as an argument for plotting with 'igraph'.
+'ggraph' natively supports the layout.
+}
+\examples{
+library(igraph)
+set.seed(12)
+g <- sample_bipartite(10, 5, "gnp", 0.5)
+fxy <- cbind(c(rep(0, 10), rep(1, 5)), NA)
+xy <- layout_with_fixed_coords(g, fxy)
+}
+\seealso{
+\link{layout_constrained_stress}
+}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 8193d3d..2d92522 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -40,6 +40,22 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// fixed_stress_major
+NumericMatrix fixed_stress_major(NumericMatrix y, NumericMatrix fixedCoords, NumericMatrix W, NumericMatrix D, int iter, double tol);
+RcppExport SEXP _graphlayouts_fixed_stress_major(SEXP ySEXP, SEXP fixedCoordsSEXP, SEXP WSEXP, SEXP DSEXP, SEXP iterSEXP, SEXP tolSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type y(ySEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type fixedCoords(fixedCoordsSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type W(WSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type D(DSEXP);
+    Rcpp::traits::input_parameter< int >::type iter(iterSEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    rcpp_result_gen = Rcpp::wrap(fixed_stress_major(y, fixedCoords, W, D, iter, tol));
+    return rcpp_result_gen;
+END_RCPP
+}
 // constrained_stress3D
 double constrained_stress3D(NumericMatrix x, NumericMatrix W, NumericMatrix D);
 RcppExport SEXP _graphlayouts_constrained_stress3D(SEXP xSEXP, SEXP WSEXP, SEXP DSEXP) {
@@ -267,6 +283,7 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_graphlayouts_constrained_stress", (DL_FUNC) &_graphlayouts_constrained_stress, 3},
     {"_graphlayouts_constrained_stress_major", (DL_FUNC) &_graphlayouts_constrained_stress_major, 6},
+    {"_graphlayouts_fixed_stress_major", (DL_FUNC) &_graphlayouts_fixed_stress_major, 6},
     {"_graphlayouts_constrained_stress3D", (DL_FUNC) &_graphlayouts_constrained_stress3D, 3},
     {"_graphlayouts_constrained_stress_major3D", (DL_FUNC) &_graphlayouts_constrained_stress_major3D, 6},
     {"_graphlayouts_criterion_angular_resolution", (DL_FUNC) &_graphlayouts_criterion_angular_resolution, 2},
diff --git a/src/constrained_stress.cpp b/src/constrained_stress.cpp
index 2613049..09b1c55 100644
--- a/src/constrained_stress.cpp
+++ b/src/constrained_stress.cpp
@@ -2,79 +2,127 @@
 using namespace Rcpp;
 
 // [[Rcpp::export]]
-double constrained_stress(NumericMatrix x, NumericMatrix W, NumericMatrix D){
-  double fct=0;
-  int n=x.nrow();
-  for(int i=0;i<(n-1);++i){
-    for(int j=(i+1);j<n;++j){
-      double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+(x(i,1)-x(j,1))*(x(i,1)-x(j,1)));
-      fct+=W(i,j)*(denom-D(i,j))*(denom-D(i,j));
+double constrained_stress(NumericMatrix x, NumericMatrix W, NumericMatrix D) {
+  double fct = 0;
+  int n = x.nrow();
+  for (int i = 0; i < (n - 1); ++i) {
+    for (int j = (i + 1); j < n; ++j) {
+      double denom = sqrt((x(i, 0) - x(j, 0)) * (x(i, 0) - x(j, 0)) +
+                          (x(i, 1) - x(j, 1)) * (x(i, 1) - x(j, 1)));
+      fct += W(i, j) * (denom - D(i, j)) * (denom - D(i, j));
     }
   }
   return fct;
 }
 
 // [[Rcpp::export]]
-NumericMatrix constrained_stress_major(NumericMatrix y,
-                           int dim,
-                           NumericMatrix W,
-                           NumericMatrix D,
-                           int iter,
-                           double tol) {
+NumericMatrix constrained_stress_major(NumericMatrix y, int dim,
+                                       NumericMatrix W, NumericMatrix D,
+                                       int iter, double tol) {
   int n = y.nrow();
+  NumericMatrix x(clone(y));
 
-  NumericMatrix x(n,2);
-  for(int i=0;i<n;++i){
-    for(int d=0;d<2;++d){
-      x(i,d)=y(i,d);
+  NumericVector wsum = rowSums(W);
+
+  double stress_old = constrained_stress(x, W, D);
+
+  for (int k = 0; k < iter; ++k) {
+    NumericMatrix xnew(n, 2);
+    xnew(_, dim - 1) = y(_, dim - 1);
+
+    for (int i = 0; i < n; ++i) {
+      for (int j = 0; j < n; ++j) {
+        if (i != j) {
+          double denom = sqrt(sum(pow(x(i, _) - x(j, _), 2)));
+          if (denom > 0.00001) {
+            int updateDim = 1 - (dim - 1);
+            xnew(i, updateDim) +=
+                W(i, j) *
+                (x(j, updateDim) +
+                 D(i, j) * (x(i, updateDim) - x(j, updateDim)) / denom);
+          }
+        }
+      }
+      int updateDim = 1 - (dim - 1);
+      xnew(i, updateDim) /= wsum[i];
     }
-  }
 
-  NumericVector wsum(n);
-  for(int i=0;i<n;++i){
-    for(int j=0;j<n;++j){
-      wsum[i]+=W(i,j);
+    double stress_new = constrained_stress(xnew, W, D);
+    double eps = (stress_old - stress_new) / stress_old;
+
+    if (eps <= tol) {
+      break;
     }
+
+    stress_old = stress_new;
+    x = clone(xnew);
   }
+  return x;
+}
 
-  double stress_old = constrained_stress(x,W,D);
+NumericMatrix replaceNA(NumericMatrix x, NumericMatrix y) {
+  int nrow = x.nrow(), ncol = x.ncol();
 
-  for(int k=0; k<iter; ++k){
-    NumericMatrix xnew(n,2); //out or in?
-    for(int i=0;i<n;++i){
-      if(dim==1){
-        xnew(i,0) = y(i,0);
-      } else{
-        xnew(i,1) = y(i,1);
+  for (int i = 0; i < nrow; ++i) {
+    for (int j = 0; j < ncol; ++j) {
+      if (!NumericMatrix::is_na(y(i, j))) {
+        x(i, j) = y(i, j);
       }
-      for(int j=0; j<n;++j){
-        if(i!=j){
-          double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+(x(i,1)-x(j,1))*(x(i,1)-x(j,1)));
-          if(denom>0.00001){
-            if(dim==2){
-              xnew(i,0) += W(i,j)*(x(j,0)+D(i,j)*(x(i,0)-x(j,0))/denom);
-            } else{
-              xnew(i,1) += W(i,j)*(x(j,1)+D(i,j)*(x(i,1)-x(j,1))/denom);
+    }
+  }
+
+  return x;
+}
+
+// [[Rcpp::export]]
+NumericMatrix fixed_stress_major(NumericMatrix y, NumericMatrix fixedCoords,
+                                 NumericMatrix W, NumericMatrix D, int iter,
+                                 double tol) {
+  int n = y.nrow();
+  NumericMatrix x(clone(y));
+
+  NumericVector wsum = rowSums(W);
+
+  double stress_old = constrained_stress(x, W, D);
+
+  for (int k = 0; k < iter; ++k) {
+    NumericMatrix xnew(n, 2);
+    std::fill(xnew.begin(), xnew.end(), 0);
+    xnew = replaceNA(xnew, fixedCoords);
+
+    for (int i = 0; i < n; ++i) {
+      for (int j = 0; j < n; ++j) {
+        if (i != j) {
+          double denom = sqrt(sum(pow(x(i, _) - x(j, _), 2)));
+          if ((denom > 0.00001)) {
+            if (!NumericMatrix::is_na(fixedCoords(i, 0))) {
+              xnew(i, 0) +=
+                  W(i, j) * (x(j, 0) + D(i, j) * (x(i, 0) - x(j, 0)) / denom);
+            }
+            if (!NumericMatrix::is_na(fixedCoords(i, 1))) {
+              xnew(i, 1) +=
+                  W(i, j) * (x(j, 1) + D(i, j) * (x(i, 1) - x(j, 1)) / denom);
             }
           }
         }
       }
-      if(dim==2){
-        xnew(i,0) = xnew(i,0)/wsum[i];
-      } else{
-        xnew(i,1) = xnew(i,1)/wsum[i];
+      if (!NumericMatrix::is_na(fixedCoords(i, 0))) {
+        xnew(i, 0) = xnew(i, 0) / wsum[i];
+      }
+      if (!NumericMatrix::is_na(fixedCoords(i, 1))) {
+        xnew(i, 1) = xnew(i, 1) / wsum[i];
       }
     }
-    double stress_new=constrained_stress(xnew,W,D);
-    double eps=(stress_old-stress_new)/stress_old;
-    if(eps<= tol){
+
+    double stress_new = constrained_stress(xnew, W, D);
+    double eps = (stress_old - stress_new) / stress_old;
+
+    if (eps <= tol) {
       break;
     }
-    stress_old=stress_new;
-    for(int i=0;i<n;++i){
-      x(i,0) = xnew(i,0);
-      x(i,1) = xnew(i,1);
-    }
+
+    stress_old = stress_new;
+    x = clone(xnew);
   }
   return x;
-}
+}
\ No newline at end of file
diff --git a/src/constrained_stress3D.cpp b/src/constrained_stress3D.cpp
index e153774..6aaaf55 100644
--- a/src/constrained_stress3D.cpp
+++ b/src/constrained_stress3D.cpp
@@ -2,103 +2,65 @@
 using namespace Rcpp;
 
 // [[Rcpp::export]]
-double constrained_stress3D(NumericMatrix x, NumericMatrix W, NumericMatrix D){
-  double fct=0;
-  int n=x.nrow();
-  for(int i=0;i<(n-1);++i){
-    for(int j=(i+1);j<n;++j){
-      double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+
-                          (x(i,1)-x(j,1))*(x(i,1)-x(j,1))+
-                          (x(i,2)-x(j,2))*(x(i,2)-x(j,2))
-      );
-      fct+=W(i,j)*(denom-D(i,j))*(denom-D(i,j));
+double constrained_stress3D(NumericMatrix x, NumericMatrix W, NumericMatrix D) {
+  double fct = 0;
+  int n = x.nrow();
+  for (int i = 0; i < (n - 1); ++i) {
+    for (int j = (i + 1); j < n; ++j) {
+      double denom = sqrt((x(i, 0) - x(j, 0)) * (x(i, 0) - x(j, 0)) +
+                          (x(i, 1) - x(j, 1)) * (x(i, 1) - x(j, 1)) +
+                          (x(i, 2) - x(j, 2)) * (x(i, 2) - x(j, 2)));
+      fct += W(i, j) * (denom - D(i, j)) * (denom - D(i, j));
     }
   }
   return fct;
 }
 
 // [[Rcpp::export]]
-NumericMatrix constrained_stress_major3D(NumericMatrix y,
-                           int dim,
-                           NumericMatrix W,
-                           NumericMatrix D,
-                           int iter,
-                           double tol) {
+NumericMatrix constrained_stress_major3D(NumericMatrix y, int dim,
+                                         NumericMatrix W, NumericMatrix D,
+                                         int iter, double tol) {
   int n = y.nrow();
+  NumericMatrix x(clone(y));
 
-  NumericMatrix x(n,3);
-  for(int i=0;i<n;++i){
-    for(int d=0;d<3;++d){
-      x(i,d)=y(i,d);
-    }
-  }
+  NumericVector wsum = rowSums(W);
 
-  NumericVector wsum(n);
-  for(int i=0;i<n;++i){
-    for(int j=0;j<n;++j){
-      wsum[i]+=W(i,j);
-    }
-  }
+  double stress_old = constrained_stress3D(x, W, D);
 
-  double stress_old = constrained_stress3D(x,W,D);
+  for (int k = 0; k < iter; ++k) {
+    NumericMatrix xnew(n, 3);
+    xnew(_, dim - 1) = y(_, dim - 1);
 
-  for(int k=0; k<iter; ++k){
-    NumericMatrix xnew(n,3); //out or in?
-    for(int i=0;i<n;++i){
-      if(dim==1){
-        xnew(i,0) = y(i,0);
-      } else if(dim==2){
-        xnew(i,1) = y(i,1);
-      } else{
-        xnew(i,2) = y(i,2);
-      }
-      for(int j=0; j<n;++j){
-        if(i!=j){
-          double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+
-                              (x(i,1)-x(j,1))*(x(i,1)-x(j,1))+
-                              (x(i,2)-x(j,2))*(x(i,2)-x(j,2))
-          );
-          if(denom>0.00001){
-            if(dim==1){
-              xnew(i,1) += W(i,j)*(x(j,1)+D(i,j)*(x(i,1)-x(j,1))/denom);
-              xnew(i,2) += W(i,j)*(x(j,2)+D(i,j)*(x(i,2)-x(j,2))/denom);
-            }
-            else if(dim==2){
-              xnew(i,0) += W(i,j)*(x(j,0)+D(i,j)*(x(i,0)-x(j,0))/denom);
-              xnew(i,2) += W(i,j)*(x(j,2)+D(i,j)*(x(i,2)-x(j,2))/denom);
-            }
-            else{
-              xnew(i,0) += W(i,j)*(x(j,0)+D(i,j)*(x(i,0)-x(j,0))/denom);
-              xnew(i,1) += W(i,j)*(x(j,1)+D(i,j)*(x(i,1)-x(j,1))/denom);
+    for (int i = 0; i < n; ++i) {
+      for (int j = 0; j < n; ++j) {
+        if (i != j) {
+          double denom = sqrt(sum(pow(x(i, _) - x(j, _), 2)));
+          if (denom > 0.00001) {
+            for (int d = 0; d < 3; ++d) {
+              if (d != dim - 1) {
+                xnew(i, d) +=
+                    W(i, j) * (x(j, d) + D(i, j) * (x(i, d) - x(j, d)) / denom);
+              }
             }
           }
         }
       }
-      if(dim==1){
-        xnew(i,1) = xnew(i,1)/wsum[i];
-        xnew(i,2) = xnew(i,2)/wsum[i];
-      }
-      else if(dim==2){
-        xnew(i,0) = xnew(i,0)/wsum[i];
-        xnew(i,2) = xnew(i,2)/wsum[i];
-      }
-      else{
-        xnew(i,0) = xnew(i,0)/wsum[i];
-        xnew(i,1) = xnew(i,1)/wsum[i];
+      for (int d = 0; d < 3; ++d) {
+        if (d != dim - 1) {
+          xnew(i, d) /= wsum[i];
+        }
       }
     }
-    double stress_new = constrained_stress3D(xnew,W,D);
-    double eps = (stress_old-stress_new)/stress_old;
 
-    if(eps<= tol){
+    double stress_new = constrained_stress3D(xnew, W, D);
+    double eps = (stress_old - stress_new) / stress_old;
+
+    if (eps <= tol) {
       break;
     }
-    stress_old=stress_new;
-    for(int i=0;i<n;++i){
-      x(i,0)=xnew(i,0);
-      x(i,1)=xnew(i,1);
-      x(i,2)=xnew(i,2);
-    }
+
+    stress_old = stress_new;
+    x = clone(xnew);
   }
   return x;
-}
+}
\ No newline at end of file