Integrating non constant vectors in curvilinear coordinates I was recently looking at some volume integrals in spherical polar coordinates in response to this question. I quickly realized that integrating using unit vectors in other coordinate systems was something I did not fully understand. So I'll ask a general question below:

DEFINITIONS:
Suppose I have some curvilinear coordinate system $(\xi_1,...,\xi_n)$in $\mathbb{R}^n$, related to the standard coordinate system by functions $\xi_k=f_k(x_1,...x_n)$. I'll denote the position vector as $$\mathbf{r} =\sum ^{n}_{k=1} x_{k}\widehat{\mathbf{e}_k}$$
And of course $\hat{\mathbf{e}_1},...\hat{\mathbf{e}_n}$ is the standard basis. Now suppose that this position vector has a different representation in curvilinear coordinates, say $$\mathbf{r} =\sum ^{n}_{k=1} g_{k}( \xi _{1} ,...,\xi _{n})\widehat{\mathbf{q}_k}$$
Where of course $$\hat{\mathbf{q}_k}=\frac{\frac{\partial \mathbf{r}}{\partial \xi _{k}}}{\| \frac{\partial \mathbf{r}}{\partial \xi _{k}} \| }$$
I'll assume they're related to the standard basis by some functions $$\hat{\mathbf{q}_k}=u_k(\hat{\mathbf{e}_1},...,\hat{\mathbf{e}_n})$$
EXAMPLE: In cylindrical coordinates this would be $\rho\hat{\boldsymbol{\rho}}+0\hat{\boldsymbol{\theta}}+z\hat{\mathbf{k}}$.

ACTUAL QUESTION:
Suppose I have some vector valued function $\mathbf{F}:\mathbb{R}^n \rightarrow \mathbb{R}^n \ ; \mathbf{F}:\mathbf{r} \mapsto \mathbf{F}(\mathbf{r}) $ and I wish to integrate this function over some region $\Omega \subseteq \mathbb{R}^n$, using curvilinear coordinates. In principle, this might look like $$\int _{\Omega }\mathbf{F}\left(\sum ^{n}_{k=1} g_{k}( \xi _{1} ,...,\xi _{n})\widehat{\mathbf{q}_k}\right)\prod ^{n}_{k=1} h_{k}\mathrm{d} \xi _{k} \\ (1)$$
Recalling the definition of our position vector $\mathbf{r}$ and using the scale factors $h_1,...,h_n$.
EXAMPLE: In polar coordinates we might want to integrate over a circle with radius $R$ centered at the origin, like this: $$\int\limits ^{2\pi }_{0}\int\limits ^{R}_{0}\mathbf{F}( r\hat{\mathbf{r}}) r\mathrm{d} r\mathrm{d} \theta \\\ (2) $$
However, here is where problems start. I can't just factor out the $\hat{\mathbf{q}_k}$'s from (1), since in general, they change depending on $\xi_1,...,\xi_n.$ The same goes for (2) - I can't just factor out the $\hat{\mathbf{r}}$ from the integral. The only good way of doing this that I know of is rewriting (2) as $$\int\limits ^{2\pi }_{0}\int\limits ^{R}_{0}\mathbf{F}( r\cos( \theta )\hat{\mathbf{i}} +r\sin( \theta )\hat{\mathbf{j}}) r\mathrm{d} r\mathrm{d} \theta $$
Or worse rewriting (1) as $$\int _{\Omega }\mathbf{F}\left(\sum ^{n}_{k=1} g_{k}( \xi _{1} ,...,\xi _{n}) u_{k}\left(\widehat{\mathbf{e}_1} ,...,\widehat{\mathbf{e}_n}\right)\right)\prod ^{n}_{k=1} h_{k}\mathrm{d} \xi _{k}$$
Then I would take care of the $\hat{\mathbf{e}_1},...,\hat{\mathbf{e}_n}$ components and then convert back into curvilinear coordinates at the end. But doesn't this kind of defeat the purpose of using the coordinate transform in the first place?!?! The whole point of what I'm doing, in principle, is integrating a function of curvilinear coordinates to return a vector in curvilinear coordinates, but it seems I have no choice but to make the integral return a vector in standard Cartesian coordinates and then go through the arduous process of converting the $\hat{\mathbf{e}_1},...,\hat{\mathbf{e}_n}$ back into $\hat{\mathbf{q}_1},...,\hat{\mathbf{q}_n}$. Am I doomed forever? Or is there some neat simplification I'm missing? Any help would be very very much appreciated.
Take care,
-K. defaoite
 A: More of a long comment than an answer:
You are correct. You are doomed. It is because you are trying to compute something that does not make sense to compute in a curvilinear coordinate (or at least I can't think of an example where it does): A vector-valued function defined as an integral over volume. 
Imagine you have the answer. How are you going to express it as in your curvilinear coordinate? The value of this function is no longer associated with a point in space (rather it is a property of the whole volume). So it does not make a lot of sense to pick a local basis set which is the property of a point in space to represent this vector in. 
Note that the position vector in your example is a very unique situation where the vector is always represented in the bases at the origin. But in most cases, when using a curvilinear coordinate, you are working with vector fields (rather than one vector space), which associate a vector space to every point in space (in more general manifolds, these are the tangent spaces of your manifold). You choose the basis of your curvilinear coordinate at that given point in space. But note that these are vectors from different vector spaces that the integral is adding up. Changing the basis as you add them up doesn't make much sense.
