Drag and lift force acting on an airfoil

When I run your code I get a FindRoot warning message:

Which makes me suspicious of the result quality. If we assume the result is correct we can speed up the integration by using the FEM for that too. We create a boundary element mesh of the foil:

bmeshFoil = 
  ToBoundaryMesh["Coordinates" -> coords[[5 ;; nn]], 
   "BoundaryElements" -> {LineElement[
      Partition[Range[Length[coords[[5 ;; nn]]]], 2, 1, 1]]}];

And integrate along the boundary:

{fdrag, flift} = 
 NIntegrate[force[{x, y}], {x, y} \[Element] bmeshFoil, 
   AccuracyGoal -> 3, PrecisionGoal -> 3] // AbsoluteTiming

(* {0.702661, {0.209457, 1.34502}} *)

Here is a partial non-NIntegrate answer that still needs work but might give you some ideas on how to proceed.

I extended the domain so that it would be easier for me to pick line segments related to the airfoil.

x1 = -2; x2 = 3; y1 = -1.5; y2 = 1.5;(*domain dimensions*)

Then I followed this example from the documentation to grab normals at line segment midpoint and the length of each segment:

bn = bmesh["BoundaryNormals"];
mean = Mean /@ GetElementCoordinates[bmesh["Coordinates"], #] & /@ 
   ElementIncidents[bmesh["BoundaryElements"]];
dist = EuclideanDistance @@@ 
     GetElementCoordinates[bmesh["Coordinates"], #] & /@ 
   ElementIncidents[bmesh["BoundaryElements"]];
ids = Flatten@
   Position[
    Flatten[mean, 1], _?(EuclideanDistance[#, {0, 0}] < 1.1 &), 1];
foilbn = bn[[1, ids]];
foilbnplt = ArrayReshape[foilbn, {1}~Join~(foilbn // Dimensions)];
foildist = dist[[1, ids]];
foildistplt = 
  ArrayReshape[foildist, {1}~Join~(foildist // Dimensions)];
foilmean = mean[[1, ids]];
foilmeanplt = 
  ArrayReshape[foilmean, {1}~Join~(foilmean // Dimensions)];
Show[bmesh["Wireframe"], 
 Graphics[MapThread[
   Arrow[{#1, #2}] &, {Join @@ foilmeanplt, 
    Join @@ (foilbnplt/5 + foilmeanplt)}]]]

It looks like we captured all the normals associated with the airfoil. You have lots of normals so I think a weighted sum should be a decent approximation to the integral.

Then, I created a function that takes a weighted sum of forces. It is fast but it needs some work and validation, but this method is similar what is done with other codes.

ClearAll[fluidStress]
fluidStress[{uif_InterpolatingFunction, vif_InterpolatingFunction, 
   pif_InterpolatingFunction}, mu_, rho_, bn_, dist_, mean_] := 
 Block[{dd, df, mesh, coords, dv, press, fx, fy, wfx, wfy, nx, ny, ux,
    uy, vx, vy}, 
  dd = Outer[(D[#1[x, y], #2]) &, {uif, vif}, {x, y}];
  df = Table[Function[{x, y}, Evaluate[dd[[i, j]]]], {i, 2}, {j, 2}];
  (*the coordinates from the foil*)
  coords = mean;
  dv = Table[df[[i, j]] @@@ coords, {i, 2}, {j, 2}];
  ux = dv[[1, 1]];
  uy = dv[[1, 2]];
  vx = dv[[2, 1]];
  vy = dv[[2, 2]];
  nx = bn[[All, 1]];
  ny = bn[[All, 2]];
  press = pif[#1, #2] & @@@ coords;
  fx = -nx*press + mu*(-2*nx*ux - ny*(uy + vx));
  fy = -ny*press + mu*(-nx*(vx + uy) - 2*ny*vy);
  wfx = dist*fx ;
  wfy = dist*fy; 
  Total /@ {wfx, wfy}
  ]
AbsoluteTiming[{fdrag, flift} = 
  fluidStress[{xVel, yVel, pressure}, 10^-3, 1, foilbn, foildist, 
   foilmean]]
(* {0.364506, {0.00244262, 0.158859}} *)