I have a rectangular domain (0<x<a, 0<y<b) and have some filled cylinders with ellipse cross-sections and different sizes which are located separately within this domain. How can I write a function (f(x,y)) describing all these ellipses in a single expression? I think I must use Heaviside function but do not know how.
Two or more ellipses can be described in a single equation by setting their products equal to zero.
$$\left(\left(\frac{x-c_{1}}{a_{1}}\right)^{2}+\left(\frac{y-d_{1}}{b_{1}}\right)^{2}-1\right)\left(\left(\frac{x-c_{2}}{a_{2}}\right)^{2}+\left(\frac{y-d_{2}}{b_{2}}\right)^{2}-1\right)=0$$
Here is a desmos graph to illustrate sucn an equation.
Desmos graph of multiple ellipses
In general $f(x,y)=0$ where
$$ f(x,y)=\prod_{k=1}^N\left(\left(\frac{x-c_{k}}{a_{k}}\right)^{2}+\left(\frac{y-d_{k}}{b_{k}}\right)^{2}-1\right) $$
ADDENDUM: OP clarifies that they want the interiors of the ellipses included as well.
I should have originally noted that what OP is asking for is an algebraic description of a subset of the coordinate plane. To couch the description in the language of functions one must specify a domain and a range. For this purpose, specify the domain as $\mathbb{R}^2$ and the range as $\{0,1\}$
Given a set $E=\{E_1,E_2,\cdots,E_N\}$ of $N$ ellipses, with equations $\left(\frac{x-c_{k}}{a_{k}}\right)^{2}+\left(\frac{y-d_{k}}{b_{k}}\right)^{2}=1$ for $0<k\le N$ define a function $f: \mathbb{R}^2\mapsto \{0,1\}$ as follows
$$ f(x,y)=\begin{cases} 1 & \text{ if } \left(\frac{x-c_{k}}{a_{k}}\right)^{2}+\left(\frac{y-d_{k}}{b_{k}}\right)^{2}\le1\text{ for some }0<k\le N\\ 0&\text{otherwise} \end{cases}$$
Then the union of the set of filled ellipses is
$$ \bigcup_{k=1}^NE_k=\{(x,y) | f(x,y)=1\} $$
RESPONSE TO OP'S MODIFIED QUESTION (dated 30 Aug 2022)
If the question is changed from 'filled ellipses' to 'filled elliptical cylinders' then it becomes a function $f: \mathbb{R}^3\mapsto \{0,1\}$ $$ f(x,y,z)=\begin{cases} 1 & \text{ if } 0\le z\le1\text{ and } \left(\frac{x-c_{k}}{a_{k}}\right)^{2}+\left(\frac{y-d_{k}}{b_{k}}\right)^{2}\le1\text{ for some }0<k\le N\\ 0&\text{otherwise} \end{cases}$$