I'm doing some buildings and, as some dimension varies, I found a way to center the doors and windows on specific faces selected via different materials. However the door is "lifted" up. More verticies on the upper edge more the door is shifted on +Z Axis. There is a way to stay at the center, the REAL center of that rectangle?
Mesh to Points—Faces node generates the points at the average of all of each face’s vertices, so when you add more vertices at the top edge of a face by dividing the neighboring face, you're adding more Z value to that position, thus moving it up. I'm not aware of a native, simpler method, so to get the real (mass) center of the face you need to do some math. Here's a setup that would calculate the center of mass for each face:
That setup's from user LordoftheFleas on BlenderArtists. Here's a slightly simplified, updated version courtesy of @Robin-Betts:
Here's the result:
I'm attaching the blend file so instead of re-creating the setup from scratch you can just append the node group "FaceCenterOfMass" and link it to your Mesh to Points node as you see in the screenshot.
Original version: 
⚠ Updated version (Blender 3.3): 
The .blend above assumes "Triangulate" node gives index of each face being triangulated to the first resulting triangle within that face. This assumption holds true for all Blender versions so far, but since version 4.0 official documentation explicitly states it is not guaranteed:
Blender 4.0 Manual: Geometry Nodes: Inspection: Geometry Randomization
The easy fix is to use the new "Sort" node in Blender 4.0, to sort triangles by captured face index. However this would make it impossible to create a tool working both in Blender 3.3 and the potential future versions changing indices. Therefore I [MvB here] decided to reimplement the setup in Blender 3.3 LTE.
2nd problem: the setup was assuming the original face area and sum of resulting tris areas are equal, but apparently it wasn't the case. Perhaps simply "Fixed" triangulation method should be used, but what I did is just accumulate the face areas the same way positions are accumulated.
Here is a second solution, that resolves pretty much all the things I was trying to do in my first solution.
Rough details of the algorithm
Like the previous solution, this solution's overall strategy is to accumulate the moment (centroid times area) and area, and to take their quotient at the end. But this solution uses a different formula for the incremental moments and areas. It does still triangulate in a sense; however it's a cheap O(n) trivial triangulation into a triangle fan from the polygon's first vertex. (Mathematically, any point in space can be used as the hub of this triangulation; if the origin is used, then the calculations amount to using the shoelace formula projected to each of the three coordinate planes. But using one of the polygon's vertices instead helps reduce numerical error in the case that the polygon is far from the origin.) The weight for each triangle's centroid is the signed area of the triangle (more precisely, the signed area of the triangle's projection onto the polygon's best-fit plane, in a sense; that is, the plane perpendicular to the polygon's area-weighted normal vector, which is a well-defined thing for any polygon, even if skew).
This formula has the following properties for various kinds of polygons.
In the case of a non-skew non-self-intersecting polygon, the answer is the usual one, mathematically agreeing with that given by the more expensive positive-triangulation-based solutions.
For a self-intersecting non-skew polygon, this definition of center-of-mass is still mathematically well-defined and triangulation independent (unlike the previous method). On the other hand, it doesn't have the nice boundedness guarantee of the previous method, and in fact it goes to infinity (in a well-defined way) as the polygon's signed area goes to zero. In the case the polygon's signed area is exactly zero (e.g. a perfect bowtie quad), the formula will yield zero divided by zero, which blender will probably turn into the zero vector (I haven't checked).
For a skew (non-coplanar) polygon, like the previous algorithm, the answer is triangulation-dependent (in this case that just means it's dependent on the choice of "first vertex" of the polygon).
More precisely: if we decompose the answer into two components, by projecting it parallel and perpendicular to the polygon's normal, then the component perpendicular to the normal is well-defined and triangulation-independent (it is, precisely, the centroid of the planar projection of the polygon) but the component parallel to the normal is triangulation-dependent, which is unfortunate. (Research problem: it would be nice to get some kind of best-fit formula for the parallel-to-normal component that is triangulation-independent even for skew polygons; I have some ideas for this.)
Here is the geometry node graph. (Thanks to Markus von Broady for simplifications.) I have added annotations wherever there is a graph edge corresponding to a concept with a clear name, so hopefully this is mostly self-documenting for anyone wanting to inspect it in detail.
The fact that this node group is now a pure field makes it extremely easy to use. For example, it can be used as drop-in replacement for a "Position" node as the Position input to a "Mesh to Points" node, without any need to thread a Geometry through it:
(Blender 4.2) 
EDIT 2024/11/26: The problem pointed out by Markus von Broady in comments (bad assumption about order of tris in Triangulation node result) has been fixed. The fix uses a safe efficient method of transfer from the triangulated geometry back to the original geometry, involving no "Nearest" transfer (bvh trees) and no "Sort Elements".
Note that this is my first answer; see my second answer for a leaner solution, which I prefer.
I started with Robin Betts' solution (given in Kuboå's answer before Markus von Broady's fix), and was able to simplify and
tune it some more.
(See Appendix at the bottom of this answer for expansion of the "Sample Inverse Index" node.)
The two "Accumulate Field" nodes are essentially the two summations in the following formula for polygon centroid (center-of-area), where the summation is over the triangles in the chosen triangulation of the polygon:
$$ \mathrm{\frac{\sum{Centroid(tri_i) Area(tri_i)}}{\sum{Area(tri_i)}}} $$
The simplifications and improvements are:
It uses only two "Capture Attribute" nodes instead of four (I'd like to use no "Capture Attribute" nodes, but I don't see how) (UPDATE: I figured it out; see my second solution)
It doesn't have a Geometry output (it seemed to be unnecessary, since it ended up being the same as the Geometry input). (I'd like to get rid of the Geometry input and all mention of Geometries, too, making this node group a "pure field", but I don't see how.) (UPDATE: I figured it out; see my second solution)
Instead of using the blender-calculated "Face Area" per-face attribute on the original faces, it calculates its own, as the sum of the triangle areas in the chosen triangulation. I think this is a good idea because it guarantees the effective per-triangle weights sum to 1, even if the polygon is slightly skew or self-intersecting (in which case this computed area will be triangulation-dependent, and in fact the thing called "Face Area" is probably computed by a different method such as the shoelace formula which will likely disagree with the sum of tri areas no matter what triangulation is chosen). That in turn guarantees the nice sanity property that the computed answer will lie within the convex hull of the input vertices, no matter how skew and snarled up the polygon is.
The screencast below demonstrates the stability improvement:
Another example of the stability improvement: when we run the improved solution on Suzanne, we can see that some of the wild answer points have been brought under control (so they are inside the convex hulls of the respective faces):
(Blender 4.2) 
This blend file has a menu that allows you to switch among various methods described in three answers and comments, including this one:
The transfer of the answer positions from the triangles back to the original faces is encapsulated in this "Sample Inverse Index" node, shown below. This is the part that both Robin Betts's solution and my solution got wrong initially.
Some notes about this part:
This method used here is safe and fast (linear-ish): its time grows proportionately to mesh size, as expected, when it is run on suzanne with catmull-clark subdivision surface with increasing subdivision levels, each of which increases the mesh size by approximately 4x. I got times L3:36ms, L4:142ms, L5:533ms, L6:2021ms.
It does not use "Nearest" sampling (referred to as the "old way" in Markus von Broady's comments, and used in the latest revision of Kuboå's answer). "Old way" seems to be shockingly slow as mesh size increases; running it (specifically MvB's latest revision of Kuboå's answer) on suzanne with subdiv surface, I got times L3:43ms, L4:371ms, L5:5131ms, L6:199343ms (more than 3 minutes, so 98.6x slower than my way at this level).
It does not use the "new way" ("Sort Elements" node, introduced fairly recently) either. I have not done any timing of "new way" yet, but my initial impression is that it's an unnecessarily heavy tool for this basic O(n) problem of inverting a lookup table.
Usage of this "Sample Inverse Index" node is analogous to the "Sample Index" node, except that type=Vector and domain=Face are hard-coded (I don't see a way to make those parameters/properties/options of a custom node), so if you want it on some other type and domain, you have to copy the node group and change these hard-codings inside it.
The next lower-level helper node is "Invert Mapping on Faces". Note that the domain=Faces is still hard-coded here, but the type=Vector is not. Expanded:
The next lower-level helper node after that is "Invert Mapping on Points". At this level, it works on a point cloud geometry, so this helper is fully reusable (no hard-coding of mesh element domain nor data type). Expanded:
Finally, the lowest-level helper is "Invert Permutation on Points", which is the core of the inverse-index-sampling problem: inverting a permutation. It works by first using a "Join Geometry" node to concatenate the input permutation with the identity permutation of the same length (using the ID attribute instead of anonymous attributes for simplicity), then applies an "Accumulate Field" node on the joined list with GroupID=Value and Value=Index (i.e. swapping the order of these two inputs, from how this node is usually used); then the second half of the "Trailing" output is the desired inverse permutation. See the annotations in the picture for more details of what's going on here.
