Grouping points by attributes from field with possibly many values

Question

I have a shapefile "cities" with one unique city_id for each city (geometry polygon, although not relevant here) :

city_name ; city_id
New-York  ; 0001
Paris     ; 0002
London    ; 0003

And a shapefile "points" (geometry single point), each with a unique point_id, and a field possible_city containing a string of at least one or more city_id (if more than one, they are separated by "-") :

point_id ; possible_city
p1       ; 0001
p2       ; 0002
p3       ; 0001-0002-0003
p4       ; 0001-0003

I am aware this is a terrible method for simulating a 1-n link compared to a true database system; however I am stuck with these shapefiles.

I would like to create an atlas for each city_id, with every point mentioning that city_id in the possible_city field. For that I need to create a coverage layer with one feature for every city_id (polygon or multipoints does not matter for me).

Following this question I tried creating a virtual layer. However I seem to be unable to link the cities.city_id and points.possible_city fields. Here is what I tried :

SELECT st_union(p.geometry) as geometry, c.city_id
FROM points p, cities c
WHERE '%'||p.possible_city||'%' LIKE c.city_id
GROUP BY c.city_id ;

I wonder if I should use a JOIN and/or a WHERE clause here. I feel this other question may be useful too, however my problem seems to lie in the fact that I am trying to link two shapefiles with two fields in this weird, non-database-approved way.

Feel free to edit the question of the post for added clarity.