Can domain decomposition methods be used without having to assemble the operator? My concern is that there is a trend towards matrix free methods for very large scale problems. I wonder how domain decomposition deal with this. I am new to DD methods, I would appreciate if anybody could point me to a paper with such implementation or if it is possible on FreeFEM++

Yes, it’s possible. You just need to define your matrix–vector product and supply a preconditioner. Then, either PETSc or HPDDM will solve the problem using your Krylov method of choice. Here, have a look at this example.