ver 1 2

Lagrangian differentiated by psi instead of rho

.gitignore

@@ -24,6 +24,9 @@ CMakeFiles/
 # Typical build directory
 /build/
 
+# VSCODE
+/.vscode/
+
 # Generated files in main directory, config/ and docs/
 /deps.mk
 config/_config.hpp
@@ -57,6 +60,11 @@ examples/ex2[0-9]
 examples/ex2[0-9]p
 examples/ex3[0-9]
 examples/ex3[0-9]p
+examples/exX
+examples/*.gf
+examples/*.mesh
+examples/*.mp4
+examples/*.txt
 
 examples/refined.mesh
 examples/displaced.mesh

data/elasticity_compliance.mesh

@@ -0,0 +1,31 @@
+MFEM mesh v1.0
+dimension
+2
+
+elements
+3
+1 3 0 1 5 4
+1 3 1 2 6 5
+1 3 2 3 7 6
+
+boundary
+8
+2 1 0 1
+1 1 4 0
+2 1 1 2
+2 1 2 3
+2 1 3 7
+2 1 5 4
+2 1 6 5
+2 1 7 6
+vertices
+8
+2
+0.000000e+00 0.000000e+00 
+1.000000e+00 0.000000e+00 
+2.000000e+00 0.000000e+00 
+3.000000e+00 0.000000e+00 
+0.000000e+00 1.000000e+00 
+1.000000e+00 1.000000e+00 
+2.000000e+00 1.000000e+00 
+3.000000e+00 1.000000e+00 

data/rect_with_top_fixed.mesh

@@ -0,0 +1,933 @@
+MFEM mesh v1.0
+dimension
+2
+
+elements
+400
+1 3 0 1 22 21
+1 3 1 2 23 22
+1 3 2 3 24 23
+1 3 3 4 25 24
+1 3 4 5 26 25
+1 3 5 6 27 26
+1 3 6 7 28 27
+1 3 7 8 29 28
+1 3 8 9 30 29
+1 3 9 10 31 30
+1 3 10 11 32 31
+1 3 11 12 33 32
+1 3 12 13 34 33
+1 3 13 14 35 34
+1 3 14 15 36 35
+1 3 15 16 37 36
+1 3 16 17 38 37
+1 3 17 18 39 38
+1 3 18 19 40 39
+1 3 19 20 41 40
+1 3 21 22 43 42
+1 3 22 23 44 43
+1 3 23 24 45 44
+1 3 24 25 46 45
+1 3 25 26 47 46
+1 3 26 27 48 47
+1 3 27 28 49 48
+1 3 28 29 50 49
+1 3 29 30 51 50
+1 3 30 31 52 51
+1 3 31 32 53 52
+1 3 32 33 54 53
+1 3 33 34 55 54
+1 3 34 35 56 55
+1 3 35 36 57 56
+1 3 36 37 58 57
+1 3 37 38 59 58
+1 3 38 39 60 59
+1 3 39 40 61 60
+1 3 40 41 62 61
+1 3 42 43 64 63
+1 3 43 44 65 64
+1 3 44 45 66 65
+1 3 45 46 67 66
+1 3 46 47 68 67
+1 3 47 48 69 68
+1 3 48 49 70 69
+1 3 49 50 71 70
+1 3 50 51 72 71
+1 3 51 52 73 72
+1 3 52 53 74 73
+1 3 53 54 75 74
+1 3 54 55 76 75
+1 3 55 56 77 76
+1 3 56 57 78 77
+1 3 57 58 79 78
+1 3 58 59 80 79
+1 3 59 60 81 80
+1 3 60 61 82 81
+1 3 61 62 83 82
+1 3 63 64 85 84
+1 3 64 65 86 85
+1 3 65 66 87 86
+1 3 66 67 88 87
+1 3 67 68 89 88
+1 3 68 69 90 89
+1 3 69 70 91 90
+1 3 70 71 92 91
+1 3 71 72 93 92
+1 3 72 73 94 93
+1 3 73 74 95 94
+1 3 74 75 96 95
+1 3 75 76 97 96
+1 3 76 77 98 97
+1 3 77 78 99 98
+1 3 78 79 100 99
+1 3 79 80 101 100
+1 3 80 81 102 101
+1 3 81 82 103 102
+1 3 82 83 104 103
+1 3 84 85 106 105
+1 3 85 86 107 106
+1 3 86 87 108 107
+1 3 87 88 109 108
+1 3 88 89 110 109
+1 3 89 90 111 110
+1 3 90 91 112 111
+1 3 91 92 113 112
+1 3 92 93 114 113
+1 3 93 94 115 114
+1 3 94 95 116 115
+1 3 95 96 117 116
+1 3 96 97 118 117
+1 3 97 98 119 118
+1 3 98 99 120 119
+1 3 99 100 121 120
+1 3 100 101 122 121
+1 3 101 102 123 122
+1 3 102 103 124 123
+1 3 103 104 125 124
+1 3 105 106 127 126
+1 3 106 107 128 127
+1 3 107 108 129 128
+1 3 108 109 130 129
+1 3 109 110 131 130
+1 3 110 111 132 131
+1 3 111 112 133 132
+1 3 112 113 134 133
+1 3 113 114 135 134
+1 3 114 115 136 135
+1 3 115 116 137 136
+1 3 116 117 138 137
+1 3 117 118 139 138
+1 3 118 119 140 139
+1 3 119 120 141 140
+1 3 120 121 142 141
+1 3 121 122 143 142
+1 3 122 123 144 143
+1 3 123 124 145 144
+1 3 124 125 146 145
+1 3 126 127 148 147
+1 3 127 128 149 148
+1 3 128 129 150 149
+1 3 129 130 151 150
+1 3 130 131 152 151
+1 3 131 132 153 152
+1 3 132 133 154 153
+1 3 133 134 155 154
+1 3 134 135 156 155
+1 3 135 136 157 156
+1 3 136 137 158 157
+1 3 137 138 159 158
+1 3 138 139 160 159
+1 3 139 140 161 160
+1 3 140 141 162 161
+1 3 141 142 163 162
+1 3 142 143 164 163
+1 3 143 144 165 164
+1 3 144 145 166 165
+1 3 145 146 167 166
+1 3 147 148 169 168
+1 3 148 149 170 169
+1 3 149 150 171 170
+1 3 150 151 172 171
+1 3 151 152 173 172
+1 3 152 153 174 173
+1 3 153 154 175 174
+1 3 154 155 176 175
+1 3 155 156 177 176
+1 3 156 157 178 177
+1 3 157 158 179 178
+1 3 158 159 180 179
+1 3 159 160 181 180
+1 3 160 161 182 181
+1 3 161 162 183 182
+1 3 162 163 184 183
+1 3 163 164 185 184
+1 3 164 165 186 185
+1 3 165 166 187 186
+1 3 166 167 188 187
+1 3 168 169 190 189
+1 3 169 170 191 190
+1 3 170 171 192 191
+1 3 171 172 193 192
+1 3 172 173 194 193
+1 3 173 174 195 194
+1 3 174 175 196 195
+1 3 175 176 197 196
+1 3 176 177 198 197
+1 3 177 178 199 198
+1 3 178 179 200 199
+1 3 179 180 201 200
+1 3 180 181 202 201
+1 3 181 182 203 202
+1 3 182 183 204 203
+1 3 183 184 205 204
+1 3 184 185 206 205
+1 3 185 186 207 206
+1 3 186 187 208 207
+1 3 187 188 209 208
+1 3 189 190 211 210
+1 3 190 191 212 211
+1 3 191 192 213 212
+1 3 192 193 214 213
+1 3 193 194 215 214
+1 3 194 195 216 215
+1 3 195 196 217 216
+1 3 196 197 218 217
+1 3 197 198 219 218
+1 3 198 199 220 219
+1 3 199 200 221 220
+1 3 200 201 222 221
+1 3 201 202 223 222
+1 3 202 203 224 223
+1 3 203 204 225 224
+1 3 204 205 226 225
+1 3 205 206 227 226
+1 3 206 207 228 227
+1 3 207 208 229 228
+1 3 208 209 230 229
+1 3 210 211 232 231
+1 3 211 212 233 232
+1 3 212 213 234 233
+1 3 213 214 235 234
+1 3 214 215 236 235
+1 3 215 216 237 236
+1 3 216 217 238 237
+1 3 217 218 239 238
+1 3 218 219 240 239
+1 3 219 220 241 240
+1 3 220 221 242 241
+1 3 221 222 243 242
+1 3 222 223 244 243
+1 3 223 224 245 244
+1 3 224 225 246 245
+1 3 225 226 247 246
+1 3 226 227 248 247
+1 3 227 228 249 248
+1 3 228 229 250 249
+1 3 229 230 251 250
+1 3 231 232 253 252
+1 3 232 233 254 253
+1 3 233 234 255 254
+1 3 234 235 256 255
+1 3 235 236 257 256
+1 3 236 237 258 257
+1 3 237 238 259 258
+1 3 238 239 260 259
+1 3 239 240 261 260
+1 3 240 241 262 261
+1 3 241 242 263 262
+1 3 242 243 264 263
+1 3 243 244 265 264
+1 3 244 245 266 265
+1 3 245 246 267 266
+1 3 246 247 268 267
+1 3 247 248 269 268
+1 3 248 249 270 269
+1 3 249 250 271 270
+1 3 250 251 272 271
+1 3 252 253 274 273
+1 3 253 254 275 274
+1 3 254 255 276 275
+1 3 255 256 277 276
+1 3 256 257 278 277
+1 3 257 258 279 278
+1 3 258 259 280 279
+1 3 259 260 281 280
+1 3 260 261 282 281
+1 3 261 262 283 282
+1 3 262 263 284 283
+1 3 263 264 285 284
+1 3 264 265 286 285
+1 3 265 266 287 286
+1 3 266 267 288 287
+1 3 267 268 289 288
+1 3 268 269 290 289
+1 3 269 270 291 290
+1 3 270 271 292 291
+1 3 271 272 293 292
+1 3 273 274 295 294
+1 3 274 275 296 295
+1 3 275 276 297 296
+1 3 276 277 298 297
+1 3 277 278 299 298
+1 3 278 279 300 299
+1 3 279 280 301 300
+1 3 280 281 302 301
+1 3 281 282 303 302
+1 3 282 283 304 303
+1 3 283 284 305 304
+1 3 284 285 306 305
+1 3 285 286 307 306
+1 3 286 287 308 307
+1 3 287 288 309 308
+1 3 288 289 310 309
+1 3 289 290 311 310
+1 3 290 291 312 311
+1 3 291 292 313 312
+1 3 292 293 314 313
+1 3 294 295 316 315
+1 3 295 296 317 316
+1 3 296 297 318 317
+1 3 297 298 319 318
+1 3 298 299 320 319
+1 3 299 300 321 320
+1 3 300 301 322 321
+1 3 301 302 323 322
+1 3 302 303 324 323
+1 3 303 304 325 324
+1 3 304 305 326 325
+1 3 305 306 327 326
+1 3 306 307 328 327
+1 3 307 308 329 328
+1 3 308 309 330 329
+1 3 309 310 331 330
+1 3 310 311 332 331
+1 3 311 312 333 332
+1 3 312 313 334 333
+1 3 313 314 335 334
+1 3 315 316 337 336
+1 3 316 317 338 337
+1 3 317 318 339 338
+1 3 318 319 340 339
+1 3 319 320 341 340
+1 3 320 321 342 341
+1 3 321 322 343 342
+1 3 322 323 344 343
+1 3 323 324 345 344
+1 3 324 325 346 345
+1 3 325 326 347 346
+1 3 326 327 348 347
+1 3 327 328 349 348
+1 3 328 329 350 349
+1 3 329 330 351 350
+1 3 330 331 352 351
+1 3 331 332 353 352
+1 3 332 333 354 353
+1 3 333 334 355 354
+1 3 334 335 356 355
+1 3 336 337 358 357
+1 3 337 338 359 358
+1 3 338 339 360 359
+1 3 339 340 361 360
+1 3 340 341 362 361
+1 3 341 342 363 362
+1 3 342 343 364 363
+1 3 343 344 365 364
+1 3 344 345 366 365
+1 3 345 346 367 366
+1 3 346 347 368 367
+1 3 347 348 369 368
+1 3 348 349 370 369
+1 3 349 350 371 370
+1 3 350 351 372 371
+1 3 351 352 373 372
+1 3 352 353 374 373
+1 3 353 354 375 374
+1 3 354 355 376 375
+1 3 355 356 377 376
+1 3 357 358 379 378
+1 3 358 359 380 379
+1 3 359 360 381 380
+1 3 360 361 382 381
+1 3 361 362 383 382
+1 3 362 363 384 383
+1 3 363 364 385 384
+1 3 364 365 386 385
+1 3 365 366 387 386
+1 3 366 367 388 387
+1 3 367 368 389 388
+1 3 368 369 390 389
+1 3 369 370 391 390
+1 3 370 371 392 391
+1 3 371 372 393 392
+1 3 372 373 394 393
+1 3 373 374 395 394
+1 3 374 375 396 395
+1 3 375 376 397 396
+1 3 376 377 398 397
+1 3 378 379 400 399
+1 3 379 380 401 400
+1 3 380 381 402 401
+1 3 381 382 403 402
+1 3 382 383 404 403
+1 3 383 384 405 404
+1 3 384 385 406 405
+1 3 385 386 407 406
+1 3 386 387 408 407
+1 3 387 388 409 408
+1 3 388 389 410 409
+1 3 389 390 411 410
+1 3 390 391 412 411
+1 3 391 392 413 412
+1 3 392 393 414 413
+1 3 393 394 415 414
+1 3 394 395 416 415
+1 3 395 396 417 416
+1 3 396 397 418 417
+1 3 397 398 419 418
+1 3 399 400 421 420
+1 3 400 401 422 421
+1 3 401 402 423 422
+1 3 402 403 424 423
+1 3 403 404 425 424
+1 3 404 405 426 425
+1 3 405 406 427 426
+1 3 406 407 428 427
+1 3 407 408 429 428
+1 3 408 409 430 429
+1 3 409 410 431 430
+1 3 410 411 432 431
+1 3 411 412 433 432
+1 3 412 413 434 433
+1 3 413 414 435 434
+1 3 414 415 436 435
+1 3 415 416 437 436
+1 3 416 417 438 437
+1 3 417 418 439 438
+1 3 418 419 440 439
+
+boundary
+80
+2 1 0 1
+2  1 21  0
+2 1 1 2
+2 1 2 3
+2 1 3 4
+2 1 4 5
+2 1 5 6
+2 1 6 7
+2 1 7 8
+2 1 8 9
+2  1  9 10
+2  1 10 11
+2  1 11 12
+2  1 12 13
+2  1 13 14
+2  1 14 15
+2  1 15 16
+2  1 16 17
+2  1 17 18
+2  1 18 19
+2  1 19 20
+2  1 20 41
+2  1 42 21
+2  1 41 62
+2  1 63 42
+2  1 62 83
+2  1 84 63
+2   1  83 104
+2   1 105  84
+2   1 104 125
+2   1 126 105
+2   1 125 146
+2   1 147 126
+2   1 146 167
+2   1 168 147
+2   1 167 188
+2   1 189 168
+2   1 188 209
+2   1 210 189
+2   1 209 230
+2   1 231 210
+2   1 230 251
+2   1 252 231
+2   1 251 272
+2   1 273 252
+2   1 272 293
+2   1 294 273
+2   1 293 314
+2   1 315 294
+2   1 314 335
+2   1 336 315
+2   1 335 356
+2   1 357 336
+2   1 356 377
+2   1 378 357
+2   1 377 398
+2   1 399 378
+2   1 398 419
+2   1 420 399
+2   1 419 440
+2   1 421 420
+2   1 422 421
+2   1 423 422
+2   1 424 423
+2   1 425 424
+2   1 426 425
+2   1 427 426
+2   1 428 427
+2   1 429 428
+1   1 430 429
+1   1 431 430
+2   1 432 431
+2   1 433 432
+2   1 434 433
+2   1 435 434
+2   1 436 435
+2   1 437 436
+2   1 438 437
+2   1 439 438
+2   1 440 439
+vertices
+441
+2
+0.000000e+00 0.000000e+00 
+1.000000e+00 0.000000e+00 
+2.000000e+00 0.000000e+00 
+3.000000e+00 0.000000e+00 
+4.000000e+00 0.000000e+00 
+5.000000e+00 0.000000e+00 
+6.000000e+00 0.000000e+00 
+7.000000e+00 0.000000e+00 
+8.000000e+00 0.000000e+00 
+9.000000e+00 0.000000e+00 
+1.000000e+01 0.000000e+00 
+1.100000e+01 0.000000e+00 
+1.200000e+01 0.000000e+00 
+1.300000e+01 0.000000e+00 
+1.400000e+01 0.000000e+00 
+1.500000e+01 0.000000e+00 
+1.600000e+01 0.000000e+00 
+1.700000e+01 0.000000e+00 
+1.800000e+01 0.000000e+00 
+1.900000e+01 0.000000e+00 
+2.000000e+01 0.000000e+00 
+0.000000e+00 1.000000e+00 
+1.000000e+00 1.000000e+00 
+2.000000e+00 1.000000e+00 
+3.000000e+00 1.000000e+00 
+4.000000e+00 1.000000e+00 
+5.000000e+00 1.000000e+00 
+6.000000e+00 1.000000e+00 
+7.000000e+00 1.000000e+00 
+8.000000e+00 1.000000e+00 
+9.000000e+00 1.000000e+00 
+1.000000e+01 1.000000e+00 
+1.100000e+01 1.000000e+00 
+1.200000e+01 1.000000e+00 
+1.300000e+01 1.000000e+00 
+1.400000e+01 1.000000e+00 
+1.500000e+01 1.000000e+00 
+1.600000e+01 1.000000e+00 
+1.700000e+01 1.000000e+00 
+1.800000e+01 1.000000e+00 
+1.900000e+01 1.000000e+00 
+2.000000e+01 1.000000e+00 
+0.000000e+00 2.000000e+00 
+1.000000e+00 2.000000e+00 
+2.000000e+00 2.000000e+00 
+3.000000e+00 2.000000e+00 
+4.000000e+00 2.000000e+00 
+5.000000e+00 2.000000e+00 
+6.000000e+00 2.000000e+00 
+7.000000e+00 2.000000e+00 
+8.000000e+00 2.000000e+00 
+9.000000e+00 2.000000e+00 
+1.000000e+01 2.000000e+00 
+1.100000e+01 2.000000e+00 
+1.200000e+01 2.000000e+00 
+1.300000e+01 2.000000e+00 
+1.400000e+01 2.000000e+00 
+1.500000e+01 2.000000e+00 
+1.600000e+01 2.000000e+00 
+1.700000e+01 2.000000e+00 
+1.800000e+01 2.000000e+00 
+1.900000e+01 2.000000e+00 
+2.000000e+01 2.000000e+00 
+0.000000e+00 3.000000e+00 
+1.000000e+00 3.000000e+00 
+2.000000e+00 3.000000e+00 
+3.000000e+00 3.000000e+00 
+4.000000e+00 3.000000e+00 
+5.000000e+00 3.000000e+00 
+6.000000e+00 3.000000e+00 
+7.000000e+00 3.000000e+00 
+8.000000e+00 3.000000e+00 
+9.000000e+00 3.000000e+00 
+1.000000e+01 3.000000e+00 
+1.100000e+01 3.000000e+00 
+1.200000e+01 3.000000e+00 
+1.300000e+01 3.000000e+00 
+1.400000e+01 3.000000e+00 
+1.500000e+01 3.000000e+00 
+1.600000e+01 3.000000e+00 
+1.700000e+01 3.000000e+00 
+1.800000e+01 3.000000e+00 
+1.900000e+01 3.000000e+00 
+2.000000e+01 3.000000e+00 
+0.000000e+00 4.000000e+00 
+1.000000e+00 4.000000e+00 
+2.000000e+00 4.000000e+00 
+3.000000e+00 4.000000e+00 
+4.000000e+00 4.000000e+00 
+5.000000e+00 4.000000e+00 
+6.000000e+00 4.000000e+00 
+7.000000e+00 4.000000e+00 
+8.000000e+00 4.000000e+00 
+9.000000e+00 4.000000e+00 
+1.000000e+01 4.000000e+00 
+1.100000e+01 4.000000e+00 
+1.200000e+01 4.000000e+00 
+1.300000e+01 4.000000e+00 
+1.400000e+01 4.000000e+00 
+1.500000e+01 4.000000e+00 
+1.600000e+01 4.000000e+00 
+1.700000e+01 4.000000e+00 
+1.800000e+01 4.000000e+00 
+1.900000e+01 4.000000e+00 
+2.000000e+01 4.000000e+00 
+0.000000e+00 5.000000e+00 
+1.000000e+00 5.000000e+00 
+2.000000e+00 5.000000e+00 
+3.000000e+00 5.000000e+00 
+4.000000e+00 5.000000e+00 
+5.000000e+00 5.000000e+00 
+6.000000e+00 5.000000e+00 
+7.000000e+00 5.000000e+00 
+8.000000e+00 5.000000e+00 
+9.000000e+00 5.000000e+00 
+1.000000e+01 5.000000e+00 
+1.100000e+01 5.000000e+00 
+1.200000e+01 5.000000e+00 
+1.300000e+01 5.000000e+00 
+1.400000e+01 5.000000e+00 
+1.500000e+01 5.000000e+00 
+1.600000e+01 5.000000e+00 
+1.700000e+01 5.000000e+00 
+1.800000e+01 5.000000e+00 
+1.900000e+01 5.000000e+00 
+2.000000e+01 5.000000e+00 
+0.000000e+00 6.000000e+00 
+1.000000e+00 6.000000e+00 
+2.000000e+00 6.000000e+00 
+3.000000e+00 6.000000e+00 
+4.000000e+00 6.000000e+00 
+5.000000e+00 6.000000e+00 
+6.000000e+00 6.000000e+00 
+7.000000e+00 6.000000e+00 
+8.000000e+00 6.000000e+00 
+9.000000e+00 6.000000e+00 
+1.000000e+01 6.000000e+00 
+1.100000e+01 6.000000e+00 
+1.200000e+01 6.000000e+00 
+1.300000e+01 6.000000e+00 
+1.400000e+01 6.000000e+00 
+1.500000e+01 6.000000e+00 
+1.600000e+01 6.000000e+00 
+1.700000e+01 6.000000e+00 
+1.800000e+01 6.000000e+00 
+1.900000e+01 6.000000e+00 
+2.000000e+01 6.000000e+00 
+0.000000e+00 7.000000e+00 
+1.000000e+00 7.000000e+00 
+2.000000e+00 7.000000e+00 
+3.000000e+00 7.000000e+00 
+4.000000e+00 7.000000e+00 
+5.000000e+00 7.000000e+00 
+6.000000e+00 7.000000e+00 
+7.000000e+00 7.000000e+00 
+8.000000e+00 7.000000e+00 
+9.000000e+00 7.000000e+00 
+1.000000e+01 7.000000e+00 
+1.100000e+01 7.000000e+00 
+1.200000e+01 7.000000e+00 
+1.300000e+01 7.000000e+00 
+1.400000e+01 7.000000e+00 
+1.500000e+01 7.000000e+00 
+1.600000e+01 7.000000e+00 
+1.700000e+01 7.000000e+00 
+1.800000e+01 7.000000e+00 
+1.900000e+01 7.000000e+00 
+2.000000e+01 7.000000e+00 
+0.000000e+00 8.000000e+00 
+1.000000e+00 8.000000e+00 
+2.000000e+00 8.000000e+00 
+3.000000e+00 8.000000e+00 
+4.000000e+00 8.000000e+00 
+5.000000e+00 8.000000e+00 
+6.000000e+00 8.000000e+00 
+7.000000e+00 8.000000e+00 
+8.000000e+00 8.000000e+00 
+9.000000e+00 8.000000e+00 
+1.000000e+01 8.000000e+00 
+1.100000e+01 8.000000e+00 
+1.200000e+01 8.000000e+00 
+1.300000e+01 8.000000e+00 
+1.400000e+01 8.000000e+00 
+1.500000e+01 8.000000e+00 
+1.600000e+01 8.000000e+00 
+1.700000e+01 8.000000e+00 
+1.800000e+01 8.000000e+00 
+1.900000e+01 8.000000e+00 
+2.000000e+01 8.000000e+00 
+0.000000e+00 9.000000e+00 
+1.000000e+00 9.000000e+00 
+2.000000e+00 9.000000e+00 
+3.000000e+00 9.000000e+00 
+4.000000e+00 9.000000e+00 
+5.000000e+00 9.000000e+00 
+6.000000e+00 9.000000e+00 
+7.000000e+00 9.000000e+00 
+8.000000e+00 9.000000e+00 
+9.000000e+00 9.000000e+00 
+1.000000e+01 9.000000e+00 
+1.100000e+01 9.000000e+00 
+1.200000e+01 9.000000e+00 
+1.300000e+01 9.000000e+00 
+1.400000e+01 9.000000e+00 
+1.500000e+01 9.000000e+00 
+1.600000e+01 9.000000e+00 
+1.700000e+01 9.000000e+00 
+1.800000e+01 9.000000e+00 
+1.900000e+01 9.000000e+00 
+2.000000e+01 9.000000e+00 
+0.000000e+00 1.000000e+01 
+1.000000e+00 1.000000e+01 
+2.000000e+00 1.000000e+01 
+3.000000e+00 1.000000e+01 
+4.000000e+00 1.000000e+01 
+5.000000e+00 1.000000e+01 
+6.000000e+00 1.000000e+01 
+7.000000e+00 1.000000e+01 
+8.000000e+00 1.000000e+01 
+9.000000e+00 1.000000e+01 
+1.000000e+01 1.000000e+01 
+1.100000e+01 1.000000e+01 
+1.200000e+01 1.000000e+01 
+1.300000e+01 1.000000e+01 
+1.400000e+01 1.000000e+01 
+1.500000e+01 1.000000e+01 
+1.600000e+01 1.000000e+01 
+1.700000e+01 1.000000e+01 
+1.800000e+01 1.000000e+01 
+1.900000e+01 1.000000e+01 
+2.000000e+01 1.000000e+01 
+0.000000e+00 1.100000e+01 
+1.000000e+00 1.100000e+01 
+2.000000e+00 1.100000e+01 
+3.000000e+00 1.100000e+01 
+4.000000e+00 1.100000e+01 
+5.000000e+00 1.100000e+01 
+6.000000e+00 1.100000e+01 
+7.000000e+00 1.100000e+01 
+8.000000e+00 1.100000e+01 
+9.000000e+00 1.100000e+01 
+1.000000e+01 1.100000e+01 
+1.100000e+01 1.100000e+01 
+1.200000e+01 1.100000e+01 
+1.300000e+01 1.100000e+01 
+1.400000e+01 1.100000e+01 
+1.500000e+01 1.100000e+01 
+1.600000e+01 1.100000e+01 
+1.700000e+01 1.100000e+01 
+1.800000e+01 1.100000e+01 
+1.900000e+01 1.100000e+01 
+2.000000e+01 1.100000e+01 
+0.000000e+00 1.200000e+01 
+1.000000e+00 1.200000e+01 
+2.000000e+00 1.200000e+01 
+3.000000e+00 1.200000e+01 
+4.000000e+00 1.200000e+01 
+5.000000e+00 1.200000e+01 
+6.000000e+00 1.200000e+01 
+7.000000e+00 1.200000e+01 
+8.000000e+00 1.200000e+01 
+9.000000e+00 1.200000e+01 
+1.000000e+01 1.200000e+01 
+1.100000e+01 1.200000e+01 
+1.200000e+01 1.200000e+01 
+1.300000e+01 1.200000e+01 
+1.400000e+01 1.200000e+01 
+1.500000e+01 1.200000e+01 
+1.600000e+01 1.200000e+01 
+1.700000e+01 1.200000e+01 
+1.800000e+01 1.200000e+01 
+1.900000e+01 1.200000e+01 
+2.000000e+01 1.200000e+01 
+0.000000e+00 1.300000e+01 
+1.000000e+00 1.300000e+01 
+2.000000e+00 1.300000e+01 
+3.000000e+00 1.300000e+01 
+4.000000e+00 1.300000e+01 
+5.000000e+00 1.300000e+01 
+6.000000e+00 1.300000e+01 
+7.000000e+00 1.300000e+01 
+8.000000e+00 1.300000e+01 
+9.000000e+00 1.300000e+01 
+1.000000e+01 1.300000e+01 
+1.100000e+01 1.300000e+01 
+1.200000e+01 1.300000e+01 
+1.300000e+01 1.300000e+01 
+1.400000e+01 1.300000e+01 
+1.500000e+01 1.300000e+01 
+1.600000e+01 1.300000e+01 
+1.700000e+01 1.300000e+01 
+1.800000e+01 1.300000e+01 
+1.900000e+01 1.300000e+01 
+2.000000e+01 1.300000e+01 
+0.000000e+00 1.400000e+01 
+1.000000e+00 1.400000e+01 
+2.000000e+00 1.400000e+01 
+3.000000e+00 1.400000e+01 
+4.000000e+00 1.400000e+01 
+5.000000e+00 1.400000e+01 
+6.000000e+00 1.400000e+01 
+7.000000e+00 1.400000e+01 
+8.000000e+00 1.400000e+01 
+9.000000e+00 1.400000e+01 
+1.000000e+01 1.400000e+01 
+1.100000e+01 1.400000e+01 
+1.200000e+01 1.400000e+01 
+1.300000e+01 1.400000e+01 
+1.400000e+01 1.400000e+01 
+1.500000e+01 1.400000e+01 
+1.600000e+01 1.400000e+01 
+1.700000e+01 1.400000e+01 
+1.800000e+01 1.400000e+01 
+1.900000e+01 1.400000e+01 
+2.000000e+01 1.400000e+01 
+0.000000e+00 1.500000e+01 
+1.000000e+00 1.500000e+01 
+2.000000e+00 1.500000e+01 
+3.000000e+00 1.500000e+01 
+4.000000e+00 1.500000e+01 
+5.000000e+00 1.500000e+01 
+6.000000e+00 1.500000e+01 
+7.000000e+00 1.500000e+01 
+8.000000e+00 1.500000e+01 
+9.000000e+00 1.500000e+01 
+1.000000e+01 1.500000e+01 
+1.100000e+01 1.500000e+01 
+1.200000e+01 1.500000e+01 
+1.300000e+01 1.500000e+01 
+1.400000e+01 1.500000e+01 
+1.500000e+01 1.500000e+01 
+1.600000e+01 1.500000e+01 
+1.700000e+01 1.500000e+01 
+1.800000e+01 1.500000e+01 
+1.900000e+01 1.500000e+01 
+2.000000e+01 1.500000e+01 
+0.000000e+00 1.600000e+01 
+1.000000e+00 1.600000e+01 
+2.000000e+00 1.600000e+01 
+3.000000e+00 1.600000e+01 
+4.000000e+00 1.600000e+01 
+5.000000e+00 1.600000e+01 
+6.000000e+00 1.600000e+01 
+7.000000e+00 1.600000e+01 
+8.000000e+00 1.600000e+01 
+9.000000e+00 1.600000e+01 
+1.000000e+01 1.600000e+01 
+1.100000e+01 1.600000e+01 
+1.200000e+01 1.600000e+01 
+1.300000e+01 1.600000e+01 
+1.400000e+01 1.600000e+01 
+1.500000e+01 1.600000e+01 
+1.600000e+01 1.600000e+01 
+1.700000e+01 1.600000e+01 
+1.800000e+01 1.600000e+01 
+1.900000e+01 1.600000e+01 
+2.000000e+01 1.600000e+01 
+0.000000e+00 1.700000e+01 
+1.000000e+00 1.700000e+01 
+2.000000e+00 1.700000e+01 
+3.000000e+00 1.700000e+01 
+4.000000e+00 1.700000e+01 
+5.000000e+00 1.700000e+01 
+6.000000e+00 1.700000e+01 
+7.000000e+00 1.700000e+01 
+8.000000e+00 1.700000e+01 
+9.000000e+00 1.700000e+01 
+1.000000e+01 1.700000e+01 
+1.100000e+01 1.700000e+01 
+1.200000e+01 1.700000e+01 
+1.300000e+01 1.700000e+01 
+1.400000e+01 1.700000e+01 
+1.500000e+01 1.700000e+01 
+1.600000e+01 1.700000e+01 
+1.700000e+01 1.700000e+01 
+1.800000e+01 1.700000e+01 
+1.900000e+01 1.700000e+01 
+2.000000e+01 1.700000e+01 
+0.000000e+00 1.800000e+01 
+1.000000e+00 1.800000e+01 
+2.000000e+00 1.800000e+01 
+3.000000e+00 1.800000e+01 
+4.000000e+00 1.800000e+01 
+5.000000e+00 1.800000e+01 
+6.000000e+00 1.800000e+01 
+7.000000e+00 1.800000e+01 
+8.000000e+00 1.800000e+01 
+9.000000e+00 1.800000e+01 
+1.000000e+01 1.800000e+01 
+1.100000e+01 1.800000e+01 
+1.200000e+01 1.800000e+01 
+1.300000e+01 1.800000e+01 
+1.400000e+01 1.800000e+01 
+1.500000e+01 1.800000e+01 
+1.600000e+01 1.800000e+01 
+1.700000e+01 1.800000e+01 
+1.800000e+01 1.800000e+01 
+1.900000e+01 1.800000e+01 
+2.000000e+01 1.800000e+01 
+0.000000e+00 1.900000e+01 
+1.000000e+00 1.900000e+01 
+2.000000e+00 1.900000e+01 
+3.000000e+00 1.900000e+01 
+4.000000e+00 1.900000e+01 
+5.000000e+00 1.900000e+01 
+6.000000e+00 1.900000e+01 
+7.000000e+00 1.900000e+01 
+8.000000e+00 1.900000e+01 
+9.000000e+00 1.900000e+01 
+1.000000e+01 1.900000e+01 
+1.100000e+01 1.900000e+01 
+1.200000e+01 1.900000e+01 
+1.300000e+01 1.900000e+01 
+1.400000e+01 1.900000e+01 
+1.500000e+01 1.900000e+01 
+1.600000e+01 1.900000e+01 
+1.700000e+01 1.900000e+01 
+1.800000e+01 1.900000e+01 
+1.900000e+01 1.900000e+01 
+2.000000e+01 1.900000e+01 
+0.000000e+00 2.000000e+01 
+1.000000e+00 2.000000e+01 
+2.000000e+00 2.000000e+01 
+3.000000e+00 2.000000e+01 
+4.000000e+00 2.000000e+01 
+5.000000e+00 2.000000e+01 
+6.000000e+00 2.000000e+01 
+7.000000e+00 2.000000e+01 
+8.000000e+00 2.000000e+01 
+9.000000e+00 2.000000e+01 
+1.000000e+01 2.000000e+01 
+1.100000e+01 2.000000e+01 
+1.200000e+01 2.000000e+01 
+1.300000e+01 2.000000e+01 
+1.400000e+01 2.000000e+01 
+1.500000e+01 2.000000e+01 
+1.600000e+01 2.000000e+01 
+1.700000e+01 2.000000e+01 
+1.800000e+01 2.000000e+01 
+1.900000e+01 2.000000e+01 
+2.000000e+01 2.000000e+01 

examples/contour_mfem_gf.m

@@ -0,0 +1,66 @@
+function H = contour_mfem_gf(ax, mesh, filename, varargin)
+args = varargin;
+% args = {};
+% mesh = read_mfem_mesh('mesh000561.mesh');
+% filename = 'rho000253.gf';
+
+dim = mesh.dim;
+fid = fopen(filename, 'r');
+str = textscan(fid, '%s', 'Delimiter', '\n');
+str = str{1};
+fes = str{2}(strfind(str{2}, ': ') + 2 : end);
+if ~startsWith(fes, 'L2')
+    error('Only L2 function is supported');
+end
+fes = fes(4:end);
+if ~startsWith(fes, 'T1')
+    error('Only Lgendre-Gauss-Lobatto basis is supported');
+end
+fes = fes(4:end);
+if ~startsWith(fes, num2str(mesh.dim))
+    error('Input solution and mesh dimension do not agree');
+end
+order = str2double(fes(end));
+data = str2double(str(6:end));
+data = reshape(data, [], mesh.nrE);
+
+if dim == 2
+    if size(mesh.v4e, 1) == 3
+        data = processData2DTri(order, data);
+        h = tricontourfh(ax, mesh, data, args{:});
+    elseif size(mesh.v4e, 1) == 4
+        data = processData2DRect(order, data);
+
+        h = tricontourfh(ax, mesh, data, args{:});
+    else
+        error('For 2D, mesh should be triangular or rectangular');
+    end
+else
+    error('Dimension must be 2D');
+end
+if nargout
+    H = h;
+end
+end
+
+function data = processData2DTri(order, data)
+
+end
+function data = processData2DRect(order, data)
+submesh = rect2rectmesh(-1, 1, -1, 1, 1, 1);
+[rtri, stri] = SNodes2D(order*2);
+r1D = SNodes1D(order);
+intpV = vertexInterp2D(rtri, stri);
+X = intpV*submesh.x4tri; X = X(:);
+Y = intpV*submesh.y4tri; Y = Y(:);
+intpBx = basisInterp1D(r1D, X, order);
+intpBy = basisInterp1D(r1D, Y, order);
+intpB = repmat(intpBx, 1, length(r1D));
+for i = 1 : length(r1D)
+    intpB(:, length(r1D)*(i-1) + 1 : length(r1D)*i) = ...
+        intpB(:, length(r1D)*(i-1) + 1 : length(r1D)*i).*intpBy(:,i);
+end
+data = intpB*data;
+data = reshape(data, length(rtri), []);
+
+end

examples/draw_solutions.m

@@ -0,0 +1,31 @@
+meshfiles = dir('mesh*.mesh');
+ufiles = dir('u*.gf');
+rhofiles = dir('rho*.gf');
+frhofiles = dir('f_rho*.gf');
+
+ax = arrayfun(@(i) axes(figure()), 1:3);
+figs = arrayfun(@(ax) ax.Parent, ax);
+recorders = arrayfun(@(name, fig) FigureRecorder(name, fig), ["u.mp4", "rho.mp4", "frho.mp4"], figs);
+for i = 1 : length(meshfiles)
+    mesh = read_mfem_mesh([meshfiles(i).folder filesep meshfiles(i).name]);
+
+    contour_mfem_gf(ax(1), mesh, [ufiles(i).folder filesep ufiles(i).name], 11);
+    title(ax(1), sprintf('Iteration %d: u', i))
+
+    contour_mfem_gf(ax(2), mesh, [rhofiles(i).folder filesep rhofiles(i).name], 0.1:0.1:0.9);
+    title(ax(2), sprintf('Iteration %d: ρ', i))
+    set(ax(2), 'Clim', [0, 1]);
+    colormap(ax(2), flipud(gray));
+
+    contour_mfem_gf(ax(3), mesh, [frhofiles(i).folder filesep frhofiles(i).name], 0.1:0.1:0.9);
+    title(ax(3), sprintf('Iteration %d: ρ̃', i))
+    set(ax(3), 'Clim', [0, 1]);
+    colormap(ax(3), flipud(gray));
+
+    arrayfun(@(a) view(a, 2), ax);
+    arrayfun(@(a) set(a, 'FontSize', 16), ax);
+    arrayfun(@(a) colorbar(a), ax);
+    arrayfun(@(a) axis(a, 'off'), ax)
+    arrayfun(@(recorder) recorder.capture(), recorders);
+end
+arrayfun(@(recorder) recorder.close(), recorders);

examples/ex36.cpp

@@ -0,0 +1,426 @@
+// Thermal compliance - Fixed Point
+//
+// min (f, u)
+// s.t -∇⋅(r(ρ̃)∇u) = f in Ω = (0, 20) × (0, 20)
+//               u = 0 on Γ = (9, 11) × {y = 20}
+//            n⋅∇u = 0 on ∂Ω \ Γ
+//        -ϵΔρ̃ + ρ̃ = ρ in Ω
+//            n⋅∇ρ̃ = 0 on ∂Ω
+//           0 ≤ ρ ≤ 1 a.e. Ω
+//
+// L = (f, u) - (r(ρ̃)∇u, ∇v) + (f, v)
+//    + (ϵ∇ρ̃, ∇λ̃) + (ρ̃, λ̃) - (S(ψ), λ)
+//    + α⁻¹D(S(ψ), S(ψ_k))
+
+#include "mfem.hpp"
+#include "proximalGalerkin.hpp"
+
+// Solution variables
+class Vars { public: enum {u, f_rho, psi, f_lam, numVars}; };
+
+void clip_abs(mfem::Vector &x, const double max_abs_val)
+{
+   for (auto &val : x) { val = std::min(max_abs_val, std::max(-max_abs_val, val)); }
+}
+
+void clip(mfem::Vector &x, const double min_val, const double max_val)
+{
+   for (auto &val : x) { val = std::min(max_val, std::max(min_val, val)); }
+}
+
+using namespace std;
+using namespace mfem;
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int problem = 0;
+   const char *mesh_file = "../data/rect_with_top_fixed.mesh";
+   int ref_levels = 2;
+   int order = 0;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   bool save_solutions = false;
+   double alpha0 = 1.0;
+   double epsilon = 1e-04;
+   double rho0 = 1e-6;
+   int simp_exp = 3;
+   double max_psi = 1e07;
+   int opt_ver = 2;
+
+   int maxit_penalty = 10000;
+   int maxit_newton = 100;
+   double tol_newton = 1e-6;
+   double tol_penalty = 1e-6;
+
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Device device(device_config);
+   device.Print();
+
+   // 2. Input data (mesh, source, ...)
+   Mesh mesh;
+   switch (problem)
+   {
+      case 0:
+         mesh = Mesh::MakeCartesian2D(32, 32, mfem::Element::QUADRILATERAL, false,
+                                      1.0, 1.0);
+         break;
+      case 1:
+         mesh = Mesh(mesh_file);
+         break;
+      default:
+         mfem_error("Undefined Problem");
+   }
+   // Mesh mesh(mesh_file);
+   int dim = mesh.Dimension();
+   const int max_attributes = mesh.bdr_attributes.Max();
+
+   double volume = 0.0;
+   for (int i=0; i<mesh.GetNE(); i++) { volume += mesh.GetElementVolume(i); }
+
+   for (int i=0; i<ref_levels; i++) { mesh.UniformRefinement(); }
+
+   // Essential boundary for each variable (numVar x numAttr)
+   Array2D<int> ess_bdr(Vars::numVars, max_attributes);
+   ess_bdr = 0;
+   switch (problem)
+   {
+      case 0:
+         ess_bdr(Vars::u, 2) = true;
+         ess_bdr(Vars::u, 3) = true;
+         break;
+      case 1:
+         ess_bdr(Vars::u, 0) = true;
+         break;
+   }
+
+   // Source and fixed temperature
+   ConstantCoefficient heat_source(1.0);
+   switch (problem)
+   {
+      case 0:
+         heat_source.constant = 1e-02;
+         break;
+      case 1:
+         heat_source.constant = 1e-03;
+         break;
+   }
+   ConstantCoefficient u_bdr(0.0);
+   const double volume_fraction = 0.4;
+   const double target_volume = volume * volume_fraction;
+
+   // 3. Finite Element Spaces and discrete solutions
+   FiniteElementSpace fes_H1_Qk2(&mesh, new H1_FECollection(order + 2, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_H1_Qk1(&mesh, new H1_FECollection(order + 1, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_H1_Qk0(&mesh, new H1_FECollection(std::max(order + 0,1),
+                                                            dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk2(&mesh, new L2_FECollection(order + 2, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk1(&mesh, new L2_FECollection(order + 1, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk0(&mesh, new L2_FECollection(order + 0, dim,
+                                                            mfem::BasisType::GaussLobatto));
+
+   Array<FiniteElementSpace*> fes(Vars::numVars);
+   fes[Vars::u] = &fes_H1_Qk1;
+   fes[Vars::f_rho] = &fes_H1_Qk1;
+   fes[Vars::psi] = &fes_L2_Qk0;
+   fes[Vars::f_lam] = fes[Vars::f_rho];
+
+   Array<int> offsets = getOffsets(fes);
+   BlockVector sol(offsets), delta_sol(offsets);
+   sol = 0.0;
+   delta_sol = 0.0;
+
+   GridFunction u(fes[Vars::u], sol.GetBlock(Vars::u));
+   GridFunction psi(fes[Vars::psi], sol.GetBlock(Vars::psi));
+   GridFunction f_rho(fes[Vars::f_rho], sol.GetBlock(Vars::f_rho));
+   GridFunction f_lam(fes[Vars::f_lam], sol.GetBlock(Vars::f_lam));
+
+   GridFunction psi_k(fes[Vars::psi]);
+
+   // Project solution
+   Array<int> ess_bdr_u;
+   ess_bdr_u.MakeRef(ess_bdr[Vars::u], max_attributes);
+   u.ProjectBdrCoefficient(u_bdr, ess_bdr_u);
+   psi = logit(volume_fraction);
+   f_rho = volume_fraction;
+
+   // 4. Define preliminary coefficients
+   ConstantCoefficient eps_cf(epsilon);
+   ConstantCoefficient alpha_k(alpha0);
+   ConstantCoefficient one_cf(1.0);
+   GridFunction zero_gf(&fes_L2_Qk2);
+   zero_gf = 0.0;
+
+   auto simp_cf = SIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto dsimp_cf = DerSIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto d2simp_cf = Der2SIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto rho_cf = SigmoidCoefficient(&psi);
+   auto rho_k_cf = SigmoidCoefficient(&psi_k);
+   auto dsigmoid_cf = DerSigmoidCoefficient(&psi);
+
+   GridFunctionCoefficient u_cf(&u);
+   GridFunctionCoefficient f_rho_cf(&f_rho);
+   GridFunctionCoefficient f_lam_cf(&f_lam);
+   GridFunctionCoefficient psi_cf(&psi);
+   GridFunctionCoefficient psi_k_cf(&psi_k);
+
+   SumCoefficient diff_rho(rho_cf, rho_k_cf, 1.0, -1.0);
+
+   GradientGridFunctionCoefficient Du(&u);
+   GradientGridFunctionCoefficient Df_rho(&f_rho);
+   GradientGridFunctionCoefficient Df_lam(&f_lam);
+
+   InnerProductCoefficient squared_normDu(Du, Du);
+
+   ProductCoefficient alph_f_lam(alpha_k, f_lam_cf);
+   ProductCoefficient alph_f_lam_dsigmoid(alph_f_lam, dsigmoid_cf);
+   ProductCoefficient psi_k_dsigmoid(psi_k_cf, dsigmoid_cf);
+   ProductCoefficient dsimp_squared_normDu(dsimp_cf, squared_normDu);
+
+   // 5. Define global system for newton iteration
+   BlockLinearSystem fixedPointSystem(offsets, fes, ess_bdr);
+   fixedPointSystem.own_blocks = true;
+   for (int i=0; i<Vars::numVars; i++)
+   {
+      fixedPointSystem.SetDiagBlockMatrix(i, new BilinearForm(fes[i]));
+   }
+
+   // Equation u
+   fixedPointSystem.GetDiagBlock(Vars::u)->AddDomainIntegrator(
+      // A += (r(ρ̃^i)∇δu, ∇v)
+      new DiffusionIntegrator(simp_cf)
+   );
+   fixedPointSystem.GetLinearForm(Vars::u)->AddDomainIntegrator(
+      new DomainLFIntegrator(heat_source)
+   );
+
+   // Equation ρ̃
+   fixedPointSystem.GetDiagBlock(Vars::f_rho)->AddDomainIntegrator(
+      new DiffusionIntegrator(eps_cf)
+   );
+   fixedPointSystem.GetDiagBlock(Vars::f_rho)->AddDomainIntegrator(
+      new MassIntegrator(one_cf)
+   );
+   fixedPointSystem.GetLinearForm(Vars::f_rho)->AddDomainIntegrator(
+      new DomainLFIntegrator(rho_cf)
+   );
+
+   switch (opt_ver)
+   {
+      case 1:
+         // Equation ψ : Version 1
+         fixedPointSystem.GetDiagBlock(Vars::psi)->AddDomainIntegrator(
+            new MassIntegrator(one_cf)
+         );
+         fixedPointSystem.GetLinearForm(Vars::psi)->AddDomainIntegrator(
+            new DomainLFIntegrator(psi_k_cf)
+         );
+         fixedPointSystem.GetLinearForm(Vars::psi)->AddDomainIntegrator(
+            new DomainLFIntegrator(alph_f_lam)
+         );
+         break;
+      case 2:
+         // Equation ψ : Version 2
+         fixedPointSystem.GetDiagBlock(Vars::psi)->AddDomainIntegrator(
+            new MassIntegrator(dsigmoid_cf)
+         );
+         fixedPointSystem.GetLinearForm(Vars::psi)->AddDomainIntegrator(
+            new DomainLFIntegrator(psi_k_dsigmoid)
+         );
+         fixedPointSystem.GetLinearForm(Vars::psi)->AddDomainIntegrator(
+            new DomainLFIntegrator(alph_f_lam_dsigmoid)
+         );
+         break;
+      default:
+         mfem_error("Undefined optimality condition version.");
+   }
+
+   // Equation λ̃
+   fixedPointSystem.GetDiagBlock(Vars::f_lam)->AddDomainIntegrator(
+      new DiffusionIntegrator(eps_cf)
+   );
+   fixedPointSystem.GetDiagBlock(Vars::f_lam)->AddDomainIntegrator(
+      new MassIntegrator(one_cf)
+   );
+   fixedPointSystem.GetLinearForm(Vars::f_lam)->AddDomainIntegrator(
+      new DomainLFIntegrator(dsimp_squared_normDu)
+   );
+
+
+
+   socketstream sout_u, sout_rho, sout_f_rho;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout_u.open(vishost, visport);
+      sout_rho.open(vishost, visport);
+      sout_f_rho.open(vishost, visport);
+      if (!sout_u)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout_u.precision(precision);
+         sout_u << "solution\n" << mesh << u;
+         sout_u << "window_title 'u'\n";
+         sout_u << "keys jmmR**c\n";
+         sout_u << flush;
+         sout_rho.precision(precision);
+         GridFunction rho(&fes_L2_Qk2);
+         rho.ProjectCoefficient(rho_cf);
+         sout_rho << "solution\n" << mesh << rho;
+         sout_rho << "window_title 'ρ'\n";
+         sout_rho << "autoscale off\n";
+         sout_rho << "valuerange 0.0 1.0\n";
+         sout_rho << "keys jmmR**c\n";
+         sout_rho << flush;
+         sout_f_rho << "solution\n" << mesh << f_rho;
+         sout_f_rho << "window_title 'ρ̃'\n";
+         sout_f_rho << "autoscale off\n";
+         sout_f_rho << "valuerange 0.0 1.0\n";
+         sout_f_rho << "keys jmmR**c\n";
+         sout_f_rho << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+
+   Array<int> ordering(0); // ordering of solving the equation
+   ordering.Append(Vars::f_rho);
+   ordering.Append(Vars::u);
+   ordering.Append(Vars::f_lam);
+   ordering.Append(Vars::psi);
+   // 6. Penalty Iteration
+   for (int k=0; k<maxit_penalty; k++)
+   {
+
+      mfem::out << "Iteration " << k + 1 << std::endl;
+      alpha_k.constant = alpha0*(k+1); // update α_k
+      psi_k = psi; // update ψ_k
+      bool newton_converged = false;
+      for (int j=0; j<maxit_newton; j++) // Newton Iteration
+      {
+         mfem::out << "\tNewton Iteration " << std::setw(5) << j + 1 << ": " <<
+                   std::flush;
+         // delta_sol = 0.0; // initialize newton difference
+         // fixedPointSystem.Assemble(delta_sol); // Update system with current solution
+         // fixedPointSystem.PCG(delta_sol); // Solve system
+         Vector old_sol(sol);
+         fixedPointSystem.Assemble(sol);
+         fixedPointSystem.PCG(sol);
+         // fixedPointSystem.SolveDiag(sol, ordering, true);
+         // Project solution
+         // NOTE: Newton stopping criteria cannot see this update. Should I consider this update?
+         const double current_volume_fraction
+            = VolumeProjection(psi, target_volume) / volume;
+         // newton successive difference
+         const double diff_newton
+            = std::sqrt(old_sol.DistanceSquaredTo(sol) / old_sol.Size());
+         mfem::out << std::scientific << diff_newton << std::endl;
+
+         if (diff_newton < tol_newton)
+         {
+            newton_converged = true;
+            break;
+         }
+      } // end of Newton iteration
+      clip_abs(psi, max_psi);
+      if (!newton_converged)
+      {
+         mfem::out << "Newton failed to converge" << std::endl;
+      }
+      if (visualization)
+      {
+
+         sout_u << "solution\n" << mesh << u << flush;
+         GridFunction rho(&fes_L2_Qk2);
+         rho.ProjectCoefficient(rho_cf);
+         sout_rho << "solution\n" << mesh << rho << "valuerange 0.0 1.0\n" << flush;
+         sout_f_rho << "solution\n" << mesh << f_rho << "valuerange 0.0 1.0\n" << flush;
+
+      }
+      if (save_solutions)
+      {
+         ostringstream filename;
+         ofstream file;
+
+         filename << "mesh" << std::setfill('0') << std::setw(6) << k << ".mesh";
+         mesh.Save(filename.str().c_str());
+         filename.str(std::string());
+
+         filename << "u" << std::setfill('0') << std::setw(6) << k << ".gf";
+         file.open(filename.str());
+         GridFunction u_high(&fes_L2_Qk2);
+         u_high.ProjectCoefficient(u_cf);
+         u_high.Save(file);
+         file.close();
+         file.clear();
+         filename.str(std::string());
+
+         GridFunction rho(&fes_L2_Qk2);
+         rho.ProjectCoefficient(rho_cf);
+         filename << "rho" << std::setfill('0') << std::setw(6) << k << ".gf";
+         file.open(filename.str());
+         rho.Save(file);
+         file.close();
+         file.clear();
+         filename.str(std::string());
+
+         GridFunction f_rho_high(&fes_L2_Qk2);
+         f_rho_high.ProjectCoefficient(f_rho_cf);
+         filename << "f_rho" << std::setfill('0') << std::setw(6) << k << ".gf";
+         file.open(filename.str());
+         f_rho_high.Save(file);
+         file.close();
+         file.clear();
+         filename.str(std::string());
+      }
+      const double diff_penalty = zero_gf.ComputeL2Error(diff_rho) /
+                                  alpha_k.constant / std::sqrt(volume);
+      mfem::out << "||ρ - ρ_k||/α_k = " << std::scientific << diff_penalty;
+      mfem::out << std::endl;
+      if (diff_penalty < tol_penalty)
+      {
+         break;
+      }
+   } // end of penalty iteration
+
+   return 0;
+}

examples/ex37.cpp

@@ -0,0 +1,482 @@
+// Thermal compliance - Newton Method
+//
+// min (f, u)
+// s.t -∇⋅(r(ρ̃)∇u) = f in Ω = (0, 20) × (0, 20)
+//               u = 0 on Γ = (9, 11) × {y = 20}
+//            n⋅∇u = 0 on ∂Ω \ Γ
+//        -ϵΔρ̃ + ρ̃ = ρ in Ω
+//            n⋅∇ρ̃ = 0 on ∂Ω
+//           0 ≤ ρ ≤ 1 a.e. Ω
+//
+// L = (f, u) - (r(ρ̃)∇u, ∇v) + (f, v)
+//    + (ϵ∇ρ̃, ∇λ̃) + (ρ̃, λ̃) - (S(ψ), λ)
+//    + α⁻¹D(S(ψ), S(ψ_k))
+
+#include "mfem.hpp"
+#include "proximalGalerkin.hpp"
+
+// Solution variables
+class Vars { public: enum {u, f_rho, psi, f_lam, numVars}; };
+
+void clip_abs(mfem::Vector &x, const double max_abs_val)
+{
+   for (auto &val : x) { val = std::min(max_abs_val, std::max(-max_abs_val, val)); }
+}
+
+void clip(mfem::Vector &x, const double min_val, const double max_val)
+{
+   for (auto &val : x) { val = std::min(max_val, std::max(min_val, val)); }
+}
+
+using namespace std;
+using namespace mfem;
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int problem = 0;
+   const char *mesh_file = "../data/rect_with_top_fixed.mesh";
+   int ref_levels = 2;
+   int order = 0;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   bool save_solutions = false;
+   double alpha0 = 1.0;
+   double epsilon = 1e-04;
+   double rho0 = 1e-6;
+   int simp_exp = 3;
+   double max_psi = 1e07;
+   int opt_ver = 2;
+
+   int maxit_penalty = 10000;
+   int maxit_newton = 100;
+   double tol_newton = 1e-6;
+   double tol_penalty = 1e-6;
+
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Device device(device_config);
+   device.Print();
+
+   // 2. Input data (mesh, source, ...)
+   Mesh mesh;
+   switch (problem)
+   {
+      case 0:
+         mesh = Mesh::MakeCartesian2D(32, 32, mfem::Element::QUADRILATERAL, false,
+                                      1.0, 1.0);
+         break;
+      case 1:
+         mesh = Mesh(mesh_file);
+         break;
+      default:
+         mfem_error("Undefined Problem");
+   }
+   // Mesh mesh(mesh_file);
+   int dim = mesh.Dimension();
+   const int max_attributes = mesh.bdr_attributes.Max();
+
+   double volume = 0.0;
+   for (int i=0; i<mesh.GetNE(); i++) { volume += mesh.GetElementVolume(i); }
+
+   for (int i=0; i<ref_levels; i++) { mesh.UniformRefinement(); }
+
+   // Essential boundary for each variable (numVar x numAttr)
+   Array2D<int> ess_bdr(Vars::numVars, max_attributes);
+   ess_bdr = 0;
+   switch (problem)
+   {
+      case 0:
+         ess_bdr(Vars::u, 2) = true;
+         ess_bdr(Vars::u, 3) = true;
+         break;
+      case 1:
+         ess_bdr(Vars::u, 0) = true;
+         break;
+   }
+
+   // Source and fixed temperature
+   ConstantCoefficient heat_source(1.0);
+   switch (problem)
+   {
+      case 0:
+         heat_source.constant = 1e-02;
+         break;
+      case 1:
+         heat_source.constant = 1e-03;
+         break;
+   }
+   ConstantCoefficient u_bdr(0.0);
+   const double volume_fraction = 0.4;
+   const double target_volume = volume * volume_fraction;
+
+   // 3. Finite Element Spaces and discrete solutions
+   FiniteElementSpace fes_H1_Qk2(&mesh, new H1_FECollection(order + 2, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_H1_Qk1(&mesh, new H1_FECollection(order + 1, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_H1_Qk0(&mesh, new H1_FECollection(std::max(order + 0,1),
+                                                            dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk2(&mesh, new L2_FECollection(order + 2, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk1(&mesh, new L2_FECollection(order + 1, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk0(&mesh, new L2_FECollection(order + 0, dim,
+                                                            mfem::BasisType::GaussLobatto));
+
+   Array<FiniteElementSpace*> fes(Vars::numVars);
+   fes[Vars::u] = &fes_H1_Qk1;
+   fes[Vars::f_rho] = &fes_H1_Qk1;
+   fes[Vars::psi] = &fes_L2_Qk0;
+   fes[Vars::f_lam] = fes[Vars::f_rho];
+
+   Array<int> offsets = getOffsets(fes);
+   BlockVector sol(offsets), delta_sol(offsets);
+   sol = 0.0;
+   delta_sol = 0.0;
+
+   GridFunction u(fes[Vars::u], sol.GetBlock(Vars::u));
+   GridFunction psi(fes[Vars::psi], sol.GetBlock(Vars::psi));
+   GridFunction f_rho(fes[Vars::f_rho], sol.GetBlock(Vars::f_rho));
+   GridFunction f_lam(fes[Vars::f_lam], sol.GetBlock(Vars::f_lam));
+
+   GridFunction psi_k(fes[Vars::psi]);
+
+   // Project solution
+   Array<int> ess_bdr_u;
+   ess_bdr_u.MakeRef(ess_bdr[Vars::u], max_attributes);
+   u.ProjectBdrCoefficient(u_bdr, ess_bdr_u);
+   psi = logit(volume_fraction);
+   f_rho = volume_fraction;
+
+   // 4. Define preliminary coefficients
+   ConstantCoefficient eps_cf(epsilon);
+   ConstantCoefficient alpha_k(alpha0);
+   ConstantCoefficient one_cf(1.0);
+   GridFunction zero_gf(&fes_L2_Qk2);
+   zero_gf = 0.0;
+
+   auto simp_cf = SIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto dsimp_cf = DerSIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto d2simp_cf = Der2SIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto rho_cf = SigmoidCoefficient(&psi);
+   auto rho_k_cf = SigmoidCoefficient(&psi_k);
+   auto dsigmoid_cf = DerSigmoidCoefficient(&psi);
+   auto d2sigmoid_cf = Der2SigmoidCoefficient(&psi);
+
+   GridFunctionCoefficient u_cf(&u);
+   GridFunctionCoefficient f_rho_cf(&f_rho);
+   GridFunctionCoefficient f_lam_cf(&f_lam);
+   GridFunctionCoefficient psi_cf(&psi);
+   GridFunctionCoefficient psi_k_cf(&psi_k);
+
+   SumCoefficient diff_rho(rho_cf, rho_k_cf, 1.0, -1.0);
+
+   GradientGridFunctionCoefficient Du(&u);
+   GradientGridFunctionCoefficient Df_rho(&f_rho);
+   GradientGridFunctionCoefficient Df_lam(&f_lam);
+
+   InnerProductCoefficient squared_normDu(Du, Du);
+
+   ProductCoefficient neg_alpha(-1.0, alpha_k);
+   ProductCoefficient alph_f_lam(alpha_k, f_lam_cf);
+   ProductCoefficient alph_f_lam_dsigmoid(alph_f_lam, dsigmoid_cf);
+   ProductCoefficient psi_k_dsigmoid(psi_k_cf, dsigmoid_cf);
+   ProductCoefficient dsimp_squared_normDu(dsimp_cf, squared_normDu);
+   ProductCoefficient d2simp_squared_normDu(d2simp_cf, squared_normDu);
+   ProductCoefficient neg_dsimp_squared_normDu(-1.0, dsimp_squared_normDu);
+   ProductCoefficient neg_d2simp_squared_normDu(-1.0, d2simp_squared_normDu);
+   ProductCoefficient neg_dsigmoid(-1.0, dsigmoid_cf);
+   ProductCoefficient neg_alpha_dsigmoid(neg_alpha, dsigmoid_cf);
+   ProductCoefficient dsimp_f_rho(dsimp_cf, f_rho_cf);
+   ProductCoefficient dsigmoid_psi_k(dsigmoid_cf, psi_k_cf);
+   ProductCoefficient neg_dsigmoid_psi(neg_dsigmoid, psi_cf);
+
+   ScalarVectorProductCoefficient dsimp_Du(dsimp_cf, Du);
+   ScalarVectorProductCoefficient dsimp_f_rho_Du(dsimp_f_rho, Du);
+
+   SumCoefficient diff_psi_and_psi_k(psi_cf, psi_k_cf, 1.0, -1.0);
+   SumCoefficient phi(diff_psi_and_psi_k, alph_f_lam, 1.0, -1.0);
+   ProductCoefficient d2sigmoid_phi(d2sigmoid_cf, phi);
+   ProductCoefficient d2sigmoid_phi_psi(d2sigmoid_phi, psi_cf);
+   SumCoefficient dsigmoid_plus_d2sigmoid_phi(dsigmoid_cf, d2sigmoid_phi);
+   ProductCoefficient neg_2_dsimp(-2.0, dsimp_cf);
+   ScalarVectorProductCoefficient neg_2_dsimp_Du(neg_2_dsimp, Du);
+   ProductCoefficient neg_d2simp_squared_normDu_f_rho(neg_d2simp_squared_normDu, f_rho_cf);
+
+   // 5. Define global system for newton iteration
+   BlockLinearSystem newtonSystem(offsets, fes, ess_bdr);
+   newtonSystem.own_blocks = true;
+   for (int i=0; i<Vars::numVars; i++)
+   {
+      newtonSystem.SetDiagBlockMatrix(i, new BilinearForm(fes[i]));
+   }
+   std::vector<std::vector<int>> offDiags
+   {
+      {Vars::u, Vars::f_rho},
+      {Vars::f_rho, Vars::psi},
+      {Vars::f_lam, Vars::u},
+      {Vars::f_lam, Vars::f_rho},
+      {Vars::psi, Vars::f_lam}
+   };
+   for (auto i : offDiags)
+   {
+      newtonSystem.SetBlockMatrix(i[0], i[1], new MixedBilinearForm(fes[i[1]], fes[i[0]]));
+   }
+
+   // Equation u
+   newtonSystem.GetDiagBlock(Vars::u)->AddDomainIntegrator(
+      // A += (r(ρ̃^i)∇δu, ∇v)
+      new DiffusionIntegrator(simp_cf)
+   );
+   newtonSystem.GetBlock(Vars::u, Vars::f_rho)->AddDomainIntegrator(
+      new TransposeIntegrator(new MixedDirectionalDerivativeIntegrator(dsimp_Du))
+   );
+   newtonSystem.GetLinearForm(Vars::u)->AddDomainIntegrator(
+      new DomainLFIntegrator(heat_source)
+   );
+   newtonSystem.GetLinearForm(Vars::u)->AddDomainIntegrator(
+      new DomainLFGradIntegrator(dsimp_f_rho_Du)
+   );
+
+   // Equation ρ̃
+   newtonSystem.GetDiagBlock(Vars::f_rho)->AddDomainIntegrator(
+      new DiffusionIntegrator(eps_cf)
+   );
+   newtonSystem.GetDiagBlock(Vars::f_rho)->AddDomainIntegrator(
+      new MassIntegrator(one_cf)
+   );
+   newtonSystem.GetBlock(Vars::f_rho, Vars::psi)->AddDomainIntegrator(
+      new MixedScalarMassIntegrator(neg_dsigmoid)
+   );
+   newtonSystem.GetLinearForm(Vars::f_rho)->AddDomainIntegrator(
+      new DomainLFIntegrator(rho_cf)
+   );
+   newtonSystem.GetLinearForm(Vars::f_rho)->AddDomainIntegrator(
+      new DomainLFIntegrator(neg_dsigmoid_psi)
+   );
+
+   switch (opt_ver)
+   {
+      case 1:
+         // Equation ψ : Version 1
+         newtonSystem.GetDiagBlock(Vars::psi)->AddDomainIntegrator(
+            new MassIntegrator(one_cf)
+         );
+         newtonSystem.GetBlock(Vars::psi, Vars::f_lam)->AddDomainIntegrator(
+            new MixedScalarMassIntegrator(neg_alpha)
+         );
+         newtonSystem.GetLinearForm(Vars::psi)->AddDomainIntegrator(
+            new DomainLFIntegrator(psi_k_cf)
+         );
+         break;
+      case 2:
+         // Equation ψ : Version 2
+         newtonSystem.GetDiagBlock(Vars::psi)->AddDomainIntegrator(
+            new MassIntegrator(dsigmoid_cf)
+         );
+         newtonSystem.GetDiagBlock(Vars::psi)->AddDomainIntegrator(
+            new MassIntegrator(d2sigmoid_phi)
+         );
+         newtonSystem.GetBlock(Vars::psi, Vars::f_lam)->AddDomainIntegrator(
+            new MixedScalarMassIntegrator(neg_alpha_dsigmoid)
+         );
+         newtonSystem.GetLinearForm(Vars::psi)->AddDomainIntegrator(
+            new DomainLFIntegrator(dsigmoid_psi_k)
+         );
+         newtonSystem.GetLinearForm(Vars::psi)->AddDomainIntegrator(
+            new DomainLFIntegrator(d2sigmoid_phi_psi)
+         );
+         break;
+      default:
+      mfem_error("Undefined optimality condition version.");
+   }
+
+   // Equation λ̃
+   newtonSystem.GetDiagBlock(Vars::f_lam)->AddDomainIntegrator(
+      new DiffusionIntegrator(eps_cf)
+   );
+   newtonSystem.GetDiagBlock(Vars::f_lam)->AddDomainIntegrator(
+      new MassIntegrator(one_cf)
+   );
+   newtonSystem.GetBlock(Vars::f_lam, Vars::f_rho)->AddDomainIntegrator(
+      new MixedScalarMassIntegrator(neg_d2simp_squared_normDu)
+   );
+   newtonSystem.GetBlock(Vars::f_lam, Vars::u)->AddDomainIntegrator(
+      new MixedDirectionalDerivativeIntegrator(neg_2_dsimp_Du)
+   );
+   newtonSystem.GetLinearForm(Vars::f_lam)->AddDomainIntegrator(
+      new DomainLFIntegrator(neg_dsimp_squared_normDu)
+   );
+   newtonSystem.GetLinearForm(Vars::f_lam)->AddDomainIntegrator(
+      new DomainLFIntegrator(neg_d2simp_squared_normDu_f_rho)
+   );
+
+
+
+   socketstream sout_u, sout_rho, sout_f_rho;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout_u.open(vishost, visport);
+      sout_rho.open(vishost, visport);
+      sout_f_rho.open(vishost, visport);
+      if (!sout_u)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout_u.precision(precision);
+         sout_u << "solution\n" << mesh << u;
+         sout_u << "window_title 'u'\n";
+         sout_u << "keys jmmR**c\n";
+         sout_u << flush;
+         sout_rho.precision(precision);
+         GridFunction rho(&fes_L2_Qk2);
+         rho.ProjectCoefficient(rho_cf);
+         sout_rho << "solution\n" << mesh << rho;
+         sout_rho << "window_title 'ρ'\n";
+         sout_rho << "autoscale off\n";
+         sout_rho << "valuerange 0.0 1.0\n";
+         sout_rho << "keys jmmR**c\n";
+         sout_rho << flush;
+         sout_f_rho << "solution\n" << mesh << f_rho;
+         sout_f_rho << "window_title 'ρ̃'\n";
+         sout_f_rho << "autoscale off\n";
+         sout_f_rho << "valuerange 0.0 1.0\n";
+         sout_f_rho << "keys jmmR**c\n";
+         sout_f_rho << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+
+   Array<int> ordering(0); // ordering of solving the equation
+   ordering.Append(Vars::f_rho);
+   ordering.Append(Vars::u);
+   ordering.Append(Vars::f_lam);
+   ordering.Append(Vars::psi);
+   // 6. Penalty Iteration
+   for (int k=0; k<maxit_penalty; k++)
+   {
+
+      mfem::out << "Iteration " << k + 1 << std::endl;
+      alpha_k.constant = alpha0*(k+1); // update α_k
+      psi_k = psi; // update ψ_k
+      bool newton_converged = false;
+      for (int j=0; j<maxit_newton; j++) // Newton Iteration
+      {
+         mfem::out << "\tNewton Iteration " << std::setw(5) << j + 1 << ": " <<
+                   std::flush;
+         Vector old_sol(sol);
+         newtonSystem.Assemble(sol);
+         newtonSystem.GMRES(sol);
+         // newton successive difference
+         const double diff_newton
+            = std::sqrt(old_sol.DistanceSquaredTo(sol) / old_sol.Size());
+         mfem::out << std::scientific << diff_newton << std::endl;
+
+         if (diff_newton < tol_newton)
+         {
+            newton_converged = true;
+            break;
+         }
+      } // end of Newton iteration
+      // Project solution
+      const double current_volume_fraction
+         = VolumeProjection(psi, target_volume) / volume;
+      clip_abs(psi, max_psi);
+      if (!newton_converged)
+      {
+         mfem::out << "Newton failed to converge" << std::endl;
+      }
+      if (visualization)
+      {
+
+         sout_u << "solution\n" << mesh << u << flush;
+         GridFunction rho(&fes_L2_Qk2);
+         rho.ProjectCoefficient(rho_cf);
+         sout_rho << "solution\n" << mesh << rho << "valuerange 0.0 1.0\n" << flush;
+         sout_f_rho << "solution\n" << mesh << f_rho << "valuerange 0.0 1.0\n" << flush;
+
+      }
+      if (save_solutions)
+      {
+         ostringstream filename;
+         ofstream file;
+
+         filename << "mesh" << std::setfill('0') << std::setw(6) << k << ".mesh";
+         mesh.Save(filename.str().c_str());
+         filename.str(std::string());
+
+         filename << "u" << std::setfill('0') << std::setw(6) << k << ".gf";
+         file.open(filename.str());
+         GridFunction u_high(&fes_L2_Qk2);
+         u_high.ProjectCoefficient(u_cf);
+         u_high.Save(file);
+         file.close();
+         file.clear();
+         filename.str(std::string());
+
+         GridFunction rho(&fes_L2_Qk2);
+         rho.ProjectCoefficient(rho_cf);
+         filename << "rho" << std::setfill('0') << std::setw(6) << k << ".gf";
+         file.open(filename.str());
+         rho.Save(file);
+         file.close();
+         file.clear();
+         filename.str(std::string());
+
+         GridFunction f_rho_high(&fes_L2_Qk2);
+         f_rho_high.ProjectCoefficient(f_rho_cf);
+         filename << "f_rho" << std::setfill('0') << std::setw(6) << k << ".gf";
+         file.open(filename.str());
+         f_rho_high.Save(file);
+         file.close();
+         file.clear();
+         filename.str(std::string());
+      }
+      const double diff_penalty = zero_gf.ComputeL2Error(diff_rho) /
+                                  alpha_k.constant / std::sqrt(volume);
+      mfem::out << "||ρ - ρ_k||/α_k = " << std::scientific << diff_penalty;
+      mfem::out << std::endl;
+      if (diff_penalty < tol_penalty)
+      {
+         break;
+      }
+   } // end of penalty iteration
+
+   return 0;
+}

examples/exX.cpp

@@ -0,0 +1,340 @@
+//                                MFEM Example X
+//
+// Compile with: make ex9
+//
+// Sample runs:
+//    exX
+//
+// Description:  This example code solves the time-dependent advection equation
+//               du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of implicit
+//               and explicit ODE time integrators, the definition of periodic
+//               boundary conditions through periodic meshes, as well as the use
+//               of GLVis for persistent visualization of a time-evolving
+//               solution. The saving of time-dependent data files for external
+//               visualization with VisIt (visit.llnl.gov) and ParaView
+//               (paraview.org) is also illustrated.
+
+#include "mfem.hpp"
+#include "proximalGalerkin.hpp"
+
+// Solution variables
+class Vars { public: enum {u, f_rho, psi, f_lam, numVars}; };
+
+void clip_abs(mfem::Vector &x, const double max_abs_val)
+{
+   for (auto &val : x) { val = std::min(max_abs_val, std::max(-max_abs_val, val)); }
+}
+
+void clip(mfem::Vector &x, const double min_val, const double max_val)
+{
+   for (auto &val : x) { val = std::min(max_val, std::max(min_val, val)); }
+}
+
+using namespace std;
+using namespace mfem;
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int problem = 0;
+   const char *mesh_file = "../data/rect_with_top_fixed.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   double alpha0 = 1.0;
+   double epsilon = 1e-03;
+   double rho0 = 1e-6;
+   int simp_exp = 3;
+   double max_psi = 1e07;
+
+   int maxit_penalty = 100;
+   int maxit_newton = 100;
+   double tol_newton = 1e-16;
+   double tol_penalty = 1e-6;
+
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Device device(device_config);
+   device.Print();
+
+   // 2. Input data (mesh, source, ...)
+   Mesh mesh(mesh_file);
+   int dim = mesh.Dimension();
+   const int max_attributes = mesh.bdr_attributes.Max();
+
+   double volume = 0.0;
+   for (int i=0; i<mesh.GetNE(); i++) { volume += mesh.GetElementVolume(i); }
+
+   for (int i=0; i<ref_levels; i++) { mesh.UniformRefinement(); }
+
+   // Essential boundary for each variable (numVar x numAttr)
+   Array2D<int> ess_bdr(Vars::numVars, max_attributes);
+   ess_bdr = 0;
+   ess_bdr(Vars::u, 0) = true;
+
+   // Source and fixed temperature
+   ConstantCoefficient heat_source(1.0);
+   ConstantCoefficient u_bdr(0.0);
+   const double volume_fraction = 0.7;
+   const double target_volume = volume * volume_fraction;
+
+   // 3. Finite Element Spaces and discrete solutions
+   FiniteElementSpace fes_H1_Qk2(&mesh, new H1_FECollection(order + 2, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_H1_Qk1(&mesh, new H1_FECollection(order + 1, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_H1_Qk0(&mesh, new H1_FECollection(order + 0, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk2(&mesh, new L2_FECollection(order + 2, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk1(&mesh, new L2_FECollection(order + 1, dim,
+                                                            mfem::BasisType::GaussLobatto));
+   FiniteElementSpace fes_L2_Qk0(&mesh, new L2_FECollection(order + 0, dim,
+                                                            mfem::BasisType::GaussLobatto));
+
+   Array<FiniteElementSpace*> fes(Vars::numVars);
+   fes[Vars::u] = &fes_H1_Qk1;
+   fes[Vars::f_rho] = &fes_H1_Qk1;
+   fes[Vars::psi] = &fes_L2_Qk0;
+   fes[Vars::f_lam] = fes[Vars::f_rho];
+
+   Array<int> offsets = getOffsets(fes);
+   BlockVector sol(offsets), delta_sol(offsets);
+   sol = 0.0;
+   delta_sol = 0.0;
+
+   GridFunction u(fes[Vars::u], sol.GetBlock(Vars::u));
+   GridFunction psi(fes[Vars::psi], sol.GetBlock(Vars::psi));
+   GridFunction f_rho(fes[Vars::f_rho], sol.GetBlock(Vars::f_rho));
+   GridFunction f_lam(fes[Vars::f_lam], sol.GetBlock(Vars::f_lam));
+
+   GridFunction psi_k(fes[Vars::psi]);
+
+   // Project solution
+   Array<int> ess_bdr_u;
+   ess_bdr_u.MakeRef(ess_bdr[Vars::u], max_attributes);
+   u.ProjectBdrCoefficient(u_bdr, ess_bdr_u);
+   psi = logit(volume_fraction);
+   f_rho = volume_fraction;
+
+   // 4. Define preliminary coefficients
+   ConstantCoefficient eps_cf(epsilon);
+   ConstantCoefficient alpha_k(alpha0);
+   ConstantCoefficient one_cf(1.0);
+
+   auto simp_cf = SIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto dsimp_cf = DerSIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto d2simp_cf = Der2SIMPCoefficient(&f_rho, simp_exp, rho0);
+   auto rho_cf = SigmoidCoefficient(&psi);
+   auto dsigmoid_cf = DerSigmoidCoefficient(&psi);
+
+   GridFunctionCoefficient u_cf(&u);
+   GridFunctionCoefficient f_rho_cf(&f_rho);
+   GridFunctionCoefficient f_lam_cf(&f_lam);
+   GridFunctionCoefficient psi_cf(&psi);
+   GridFunctionCoefficient psi_k_cf(&psi_k);
+
+   GradientGridFunctionCoefficient Du(&u);
+   GradientGridFunctionCoefficient Df_rho(&f_rho);
+   GradientGridFunctionCoefficient Df_lam(&f_lam);
+
+   InnerProductCoefficient squared_normDu(Du, Du);
+
+   ProductCoefficient alph_f_lam(alpha_k, f_lam_cf);
+   ProductCoefficient neg_f_lam(-1.0, f_lam_cf);
+   ProductCoefficient neg_simp(-1.0, simp_cf);
+   ProductCoefficient neg_dsimp(-1.0, dsimp_cf);
+   ProductCoefficient dsimp_times2(2.0, dsimp_cf);
+   ProductCoefficient neg_dsimp_squared_normDu(neg_dsimp, squared_normDu);
+   ProductCoefficient d2simp_squared_normDu(d2simp_cf, squared_normDu);
+   ProductCoefficient neg_dsigmoid(-1.0, dsigmoid_cf);
+
+   ScalarVectorProductCoefficient neg_simp_Du(neg_simp, Du);
+   ScalarVectorProductCoefficient neg_eps_Df_rho(-epsilon, Df_rho);
+   ScalarVectorProductCoefficient neg_eps_Df_lam(-epsilon, Df_lam);
+   ScalarVectorProductCoefficient dsimp_Du(dsimp_cf, Du);
+   ScalarVectorProductCoefficient dsimp_Du_times2(dsimp_times2, Du);
+
+   SumCoefficient diff_filter(rho_cf, f_rho_cf, 1.0, -1.0);
+   SumCoefficient diff_psi_k(psi_k_cf, psi_cf, 1.0, -1.0);
+   SumCoefficient diff_psi_grad(diff_psi_k, alph_f_lam, 1.0, -1.0);
+
+   // 5. Define global system for newton iteration
+   BlockLinearSystem newtonSystem(offsets, fes, ess_bdr);
+   newtonSystem.own_blocks = true;
+   for (int i=0; i<Vars::numVars; i++)
+   {
+      newtonSystem.SetDiagBlockMatrix(i, new BilinearForm(fes[i]));
+   }
+   std::vector<std::vector<int>> offDiagBlocks
+   {
+      {Vars::u, Vars::f_rho},
+      {Vars::f_rho, Vars::psi},
+      {Vars::psi, Vars::f_lam},
+      {Vars::f_lam, Vars::u},
+      {Vars::f_lam, Vars::f_rho}
+   };
+   for (auto idx: offDiagBlocks)
+   {
+      newtonSystem.SetBlockMatrix(idx[0], idx[1], new MixedBilinearForm(fes[idx[1]],
+                                                                        fes[idx[0]]));
+   }
+
+
+   // Equation u
+   newtonSystem.GetDiagBlock(Vars::u)->AddDomainIntegrator(
+      // A += (r(ρ̃^i)∇δu, ∇v)
+      new DiffusionIntegrator(simp_cf)
+   );
+   newtonSystem.GetBlock(Vars::u, Vars::f_rho)->AddDomainIntegrator(
+      // A += ((r'(ρ̃^i)∇u) δρ̃, ∇v)
+      new TransposeIntegrator(new MixedDirectionalDerivativeIntegrator(dsimp_Du))
+   );
+   newtonSystem.GetLinearForm(Vars::u)->AddDomainIntegrator(
+      // b += (f, v)
+      new DomainLFIntegrator(heat_source)
+   );
+   newtonSystem.GetLinearForm(Vars::u)->AddDomainIntegrator(
+      // b += -(r(ρ̃^i)∇u^i, ∇v)
+      new DomainLFGradIntegrator(neg_simp_Du)
+   );
+
+   // Equation ρ̃
+   newtonSystem.GetDiagBlock(Vars::f_rho)->AddDomainIntegrator(
+      // A += (ϵ∇δρ̃, ∇μ̃)
+      new DiffusionIntegrator(eps_cf)
+   );
+   newtonSystem.GetDiagBlock(Vars::f_rho)->AddDomainIntegrator(
+      // A += (δρ̃, μ̃)
+      new MassIntegrator()
+   );
+   newtonSystem.GetBlock(Vars::f_rho, Vars::psi)->AddDomainIntegrator(
+      // A += -(sig'(ψ^i)δψ, μ̃)
+      new MixedScalarMassIntegrator(neg_dsigmoid)
+   );
+   newtonSystem.GetLinearForm(Vars::f_rho)->AddDomainIntegrator(
+      // b += -(ϵ∇ρ̃^i, ∇μ̃)
+      new DomainLFGradIntegrator(neg_eps_Df_rho)
+   );
+   newtonSystem.GetLinearForm(Vars::f_rho)->AddDomainIntegrator(
+      // b += (ρ^i-ρ̃^i, μ̃)
+      new DomainLFIntegrator(diff_filter)
+   );
+
+   // Equation ψ
+   newtonSystem.GetDiagBlock(Vars::psi)->AddDomainIntegrator(
+      // A += (δψ, φ)
+      new MassIntegrator()
+   );
+   newtonSystem.GetBlock(Vars::psi, Vars::f_lam)->AddDomainIntegrator(
+      // A += (α_k δλ̃, φ)
+      new MixedScalarMassIntegrator(alpha_k)
+   );
+   newtonSystem.GetLinearForm(Vars::psi)->AddDomainIntegrator(
+      // b += (ψ_k - ψ^i - α_k λ̃^i, φ)
+      new DomainLFIntegrator(diff_psi_grad)
+   );
+
+   // Equation f_lam
+   newtonSystem.GetDiagBlock(Vars::f_lam)->AddDomainIntegrator(
+      // A += (ϵ∇δλ̃, μ̃)
+      new DiffusionIntegrator(eps_cf)
+   );
+   newtonSystem.GetDiagBlock(Vars::f_lam)->AddDomainIntegrator(
+      // A += (δλ̃, μ̃)
+      new MassIntegrator()
+   );
+   newtonSystem.GetBlock(Vars::f_lam, Vars::u)->AddDomainIntegrator(
+      // A += (2r'(ρ̃)∇δu, μ̃)
+      new MixedDirectionalDerivativeIntegrator(dsimp_Du_times2)
+   );
+   newtonSystem.GetBlock(Vars::f_lam, Vars::f_rho)->AddDomainIntegrator(
+      // A += (r''(ρ̃)||∇u||^2 δρ̃, μ̃)
+      new MixedScalarMassIntegrator(d2simp_squared_normDu)
+   );
+   newtonSystem.GetLinearForm(Vars::f_lam)->AddDomainIntegrator(
+      // b += -(r'(ρ̃^i)||∇u^i||^2, μ̃)
+      new DomainLFIntegrator(neg_dsimp_squared_normDu)
+   );
+   newtonSystem.GetLinearForm(Vars::f_lam)->AddDomainIntegrator(
+      // b += -(ϵ∇λ̃^i, ∇μ̃)
+      new DomainLFGradIntegrator(neg_eps_Df_lam)
+   );
+   newtonSystem.GetLinearForm(Vars::f_lam)->AddDomainIntegrator(
+      // b += -(λ̃, μ̃)
+      new DomainLFIntegrator(neg_f_lam)
+   );
+
+   // 6. Penalty Iteration
+   for (int k=0; k<maxit_penalty; k++)
+   {
+      mfem::out << "Iteration " << k + 1 << std::endl;
+      alpha_k.constant = alpha0*(k+1); // update α_k
+      psi_k = psi; // update ψ_k
+      for (int j=0; j<maxit_newton; j++) // Newton Iteration
+      {
+         mfem::out << "\tNewton Iteration " << j + 1 << ": ";
+         delta_sol = 0.0; // initialize newton difference
+         newtonSystem.Assemble(delta_sol); // Update system with current solution
+         newtonSystem.GMRES(delta_sol); // Solve system
+         sol += delta_sol; // Update solution
+         // newton successive difference
+         const double diff_newton = std::sqrt(
+                                       std::pow(delta_sol.Norml2(), 2) / delta_sol.Size()
+                                    );
+         // Project solution
+         // NOTE: Newton stopping criteria cannot see this update. Should I consider this update?
+         const double current_volume_fraction = VolumeProjection(psi,
+                                                                 target_volume) / volume;
+         clip_abs(psi, max_psi);
+         mfem::out << std::scientific << diff_newton << ", ∫ρ / |Ω| = " << std::fixed
+                   << current_volume_fraction << std::endl;
+
+         if (diff_newton < tol_newton)
+         {
+            break;
+         }
+      } // end of Newton iteration
+      const double diff_penalty = std::sqrt(
+                                     psi_k.DistanceSquaredTo(psi) / delta_sol.Size()
+                                  );
+      mfem::out << "||ψ - ψ_k|| = " << std::scientific << diff_penalty << std::endl;
+      if (diff_penalty < tol_penalty)
+      {
+         break;
+      }
+   }
+
+   return 0;
+}
+

examples/makefile

@@ -23,7 +23,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
   ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-  ex31 ex33
+  ex31 ex33 exX
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
   ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
   ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p

examples/plot_mfem_gf.m

@@ -0,0 +1,86 @@
+function H = plot_mfem_gf(ax, mesh, filename, varargin)
+args = varargin;
+% args = {};
+% mesh = read_mfem_mesh('mesh000561.mesh');
+% filename = 'rho000253.gf';
+
+dim = mesh.dim;
+fid = fopen(filename, 'r');
+str = textscan(fid, '%s', 'Delimiter', '\n');
+str = str{1};
+fes = str{2}(strfind(str{2}, ': ') + 2 : end);
+if ~startsWith(fes, 'L2')
+    error('Only L2 function is supported');
+end
+fes = fes(4:end);
+if ~startsWith(fes, 'T1')
+    error('Only Lgendre-Gauss-Lobatto basis is supported');
+end
+fes = fes(4:end);
+if ~startsWith(fes, num2str(mesh.dim))
+    error('Input solution and mesh dimension do not agree');
+end
+order = str2double(fes(end));
+data = str2double(str(6:end));
+data = reshape(data, [], mesh.nrE);
+
+if dim == 1
+    data = processData1D(order, data);
+    r1D = SNodes1D(order);
+    h = plotlines(ax, mesh.x4tri*r1D, data, args{:});
+elseif dim == 2
+    if size(mesh.v4e, 1) == 3
+        data = processData2DTri(order, data);
+        h = trisurfh(ax, mesh, data, args{:});
+    elseif size(mesh.v4e, 1) == 4
+        data = processData2DRect(order, data);
+
+        h = polysurfh(ax, mesh, data, args{:});
+        if length(h) > 1
+            delete(h(end));
+        end
+    else
+        error('For 2D, mesh should be triangular or rectangular');
+    end
+elseif dim == 3
+    error('3D visualization not supported yet.');
+%     if size(mesh.v4e, 1) == 4
+%         data = processData3DTet(order, data);
+%     elseif size(mesh.v4e, 1) == 8
+%         data = processData3DCube(order, data);
+%     else
+%         error('For 3D, mesh should be tetrahedral or cubical');
+%     end
+else
+    error('Dimension is higher than 3');
+end
+if nargout
+    H = h;
+end
+end
+
+function data = processData1D(order, data)
+
+end
+
+function data = processData2DTri(order, data)
+
+end
+function data = processData2DRect(order, data)
+submesh = rect2rectmesh(-1, 1, -1, 1, 1, 1);
+[rtri, stri] = SNodes2D(order*2);
+r1D = SNodes1D(order);
+intpV = vertexInterp2D(rtri, stri);
+X = intpV*submesh.x4tri; X = X(:);
+Y = intpV*submesh.y4tri; Y = Y(:);
+intpBx = basisInterp1D(r1D, X, order);
+intpBy = basisInterp1D(r1D, Y, order);
+intpB = repmat(intpBx, 1, length(r1D));
+for i = 1 : length(r1D)
+    intpB(:, length(r1D)*(i-1) + 1 : length(r1D)*i) = ...
+        intpB(:, length(r1D)*(i-1) + 1 : length(r1D)*i).*intpBy(:,i);
+end
+data = intpB*data;
+data = reshape(data, length(rtri), []);
+
+end

examples/proximalGalerkin.hpp

@@ -0,0 +1,591 @@
+#ifndef PROXIMAL_GALERKIN_HPP
+#define PROXIMAL_GALERKIN_HPP
+#include "mfem.hpp"
+
+
+double sigmoid(const double x)
+{
+   if (x < 0)
+   {
+      const double exp_x = std::exp(x);
+      return exp_x / (1.0 + exp_x);
+   }
+   return 1.0 / (1.0 + std::exp(-x));
+}
+
+double dsigmoiddx(const double x)
+{
+   const double tmp = sigmoid(x);
+   return tmp*(1 - tmp);
+}
+
+double d2sigmoiddx2(const double x)
+{
+   const double tmp = sigmoid(x);
+   return tmp*(1 - tmp)*(1-2*tmp);
+}
+
+double logit(const double x)
+{
+   return std::log(x / (1.0 - x));
+}
+
+double simpRule(const double rho, const int exponent, const double rho0)
+{
+   return rho0 + (1.0 - rho0)*std::pow(rho, exponent);
+}
+
+double dsimpRuledx(const double rho, const int exponent, const double rho0)
+{
+   return exponent*(1.0 - rho0)*std::pow(rho, exponent - 1);
+}
+
+double d2simpRuledx2(const double rho, const int exponent, const double rho0)
+{
+   return exponent*(exponent - 1)*(1.0 - rho0)*std::pow(rho, exponent - 2);
+}
+
+
+namespace mfem
+{
+class MultiProductCoefficient : public Coefficient
+{
+private:
+   double aConst;
+   Array<Coefficient *> a;
+
+public:
+   /// Constructor with one coefficient.  Result is A * B.
+   MultiProductCoefficient(double A, Coefficient &B)
+      : aConst(A), a(1) { a[0] = &B;}
+
+   /// Constructor with two coefficients.  Result is A * B.
+   MultiProductCoefficient(Coefficient &A, Coefficient &B)
+      : aConst(1.0), a(2) { a[0] = &A; a[1] = &B;}
+
+   /// Constructor with two coefficients.  Result is A * B.
+   MultiProductCoefficient(double A, Coefficient &B, Coefficient &C)
+      : aConst(A), a(2) { a[0] = &B; a[1] = &C;}
+
+   /// Constructor with two coefficients.  Result is A * B.
+   MultiProductCoefficient(Coefficient &A, Coefficient &B, Coefficient &C)
+      : aConst(1.0), a(3) { a[0] = &A; a[1] = &B; a[2] = &C; }
+
+   /// Constructor with two coefficients.  Result is A * B.
+   MultiProductCoefficient(double A, Coefficient &B, Coefficient &C,
+                           Coefficient &D)
+      : aConst(A), a(3) { a[0] = &B; a[1] = &C; a[2] = &D;}
+
+   /// Constructor with two coefficients.  Result is A * B.
+   MultiProductCoefficient(Array<Coefficient *> &A)
+      : aConst(1.0), a(A) {}
+   /// Constructor with two coefficients.  Result is A * B.
+   MultiProductCoefficient(double A, Array<Coefficient *> &B)
+      : aConst(A), a(B) {}
+
+   /// Set the time for internally stored coefficients
+   void SetTime(double t)
+   {
+      for (auto &c : a) { c->SetTime(t); }
+   }
+
+   void Mult(double A) { aConst *= A; }
+   void Mult(Coefficient &A) { a.Append(&A); }
+
+   /// Evaluate the coefficient at @a ip.
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      double val = aConst;
+      for (auto &c : a) { val *= c->Eval(T, ip); }
+      return val;
+   }
+};
+class MultiProductVectorCoefficient : public VectorCoefficient
+{
+private:
+   double aConst;
+   Array<Coefficient *> a;
+   VectorCoefficient * v;
+
+public:
+   /// Constructor with one coefficient.  Result is A * B.
+   MultiProductVectorCoefficient(double A, VectorCoefficient &V)
+      : aConst(A), a(0), v(&V), VectorCoefficient(V.GetVDim()) {}
+
+   MultiProductVectorCoefficient(Coefficient &A, VectorCoefficient &V)
+      : aConst(1.0), a(0), v(&V), VectorCoefficient(V.GetVDim()) { a.Append(&A); }
+   /// Constructor with one coefficient.  Result is A * V.
+   MultiProductVectorCoefficient(double A, Coefficient &B, VectorCoefficient &V)
+      : aConst(A), a(0), v(&V), VectorCoefficient(V.GetVDim()) { a.Append(&B); }
+
+   MultiProductVectorCoefficient(Coefficient &A, Coefficient &B,
+                                 VectorCoefficient &V)
+      : aConst(1.0), a(0), v(&V), VectorCoefficient(V.GetVDim()) { a.Append(&A); a.Append(&B); }
+
+   /// Constructor with one coefficient.  Result is A * V.
+   MultiProductVectorCoefficient(double A, Coefficient &B, Coefficient &C,
+                                 VectorCoefficient &V)
+      : aConst(A), a(0), v(&V), VectorCoefficient(V.GetVDim()) { a.Append(&B); a.Append(&C); }
+
+   MultiProductVectorCoefficient(Coefficient &A, Coefficient &B, Coefficient &C,
+                                 VectorCoefficient &V)
+      : aConst(1.0), a(0), v(&V), VectorCoefficient(V.GetVDim()) { a.Append(&A); a.Append(&B); a.Append(&C); }
+
+   MultiProductVectorCoefficient(Array<Coefficient*> &A, VectorCoefficient &V)
+      : aConst(1.0), a(A), v(&V), VectorCoefficient(V.GetVDim()) {}
+
+   MultiProductVectorCoefficient(double A, Array<Coefficient*> &B,
+                                 VectorCoefficient &V)
+      : aConst(A), a(B), v(&V), VectorCoefficient(V.GetVDim()) {}
+
+   /// Set the time for internally stored coefficients
+   void SetTime(double t)
+   {
+      for (auto &c : a) { c->SetTime(t); }
+      v->SetTime(t);
+   };
+
+   void Mult(double A) { aConst *= A; }
+   void Mult(Coefficient &A) { a.Append(&A); }
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      double val = aConst;
+      for (auto &c : a) { val *= c->Eval(T, ip); }
+      v->Eval(V, T, ip);
+      V *= val;
+   };
+   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                     const IntegrationRule &ir)
+   {
+      const int N = ir.GetNPoints();
+      M.SetSize(v->GetVDim(), N);
+      Vector col;
+      for (int i=0; i<N; i++)
+      {
+         M.GetColumnReference(i, col);
+         const IntegrationPoint &ip = ir.IntPoint(i);
+         T.SetIntPoint(&ip);
+         v->Eval(col, T, ip);
+         double val = aConst;
+         for (auto &c : a) { val *= c->Eval(T, ip); }
+         col *= val;
+      }
+   }
+};
+
+
+class MappedGridFunctionCoefficient :public GridFunctionCoefficient
+{
+   typedef std::__1::function<double (const double)> Mapping;
+public:
+   MappedGridFunctionCoefficient(GridFunction *gf, Mapping f_x, int comp=1)
+      :GridFunctionCoefficient(gf, comp), map(f_x) { }
+
+   /// Evaluate the coefficient at @a ip.
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      return map(GridFunctionCoefficient::Eval(T, ip));
+   }
+
+   inline void SetMapping(Mapping &f_x) { map = f_x; }
+private:
+   Mapping map;
+};
+
+MappedGridFunctionCoefficient SIMPCoefficient(GridFunction *gf,
+                                              const int exponent, const double rho0)
+{
+   auto map = [exponent, rho0](const double x) {return simpRule(x, exponent, rho0); };
+   return MappedGridFunctionCoefficient(gf, map);
+}
+
+MappedGridFunctionCoefficient DerSIMPCoefficient(GridFunction *gf,
+                                                 const int exponent, const double rho0)
+{
+   auto map = [exponent, rho0](const double x) {return dsimpRuledx(x, exponent, rho0); };
+   return MappedGridFunctionCoefficient(gf, map);
+}
+
+MappedGridFunctionCoefficient Der2SIMPCoefficient(GridFunction *gf,
+                                                  const int exponent, const double rho0)
+{
+   auto map = [exponent, rho0](const double x) {return d2simpRuledx2(x, exponent, rho0); };
+   return MappedGridFunctionCoefficient(gf, map);
+}
+
+MappedGridFunctionCoefficient SigmoidCoefficient(GridFunction *gf)
+{
+   return MappedGridFunctionCoefficient(gf, sigmoid);
+}
+
+MappedGridFunctionCoefficient DerSigmoidCoefficient(GridFunction *gf)
+{
+   return MappedGridFunctionCoefficient(gf, dsigmoiddx);
+}
+
+MappedGridFunctionCoefficient Der2SigmoidCoefficient(GridFunction *gf)
+{
+   return MappedGridFunctionCoefficient(gf, d2sigmoiddx2);
+}
+
+double VolumeProjection(GridFunction &psi, const double target_volume,
+                        const double tol=1e-12,
+                        const int max_its=10)
+{
+   auto sigmoid_psi = SigmoidCoefficient(&psi);
+   auto der_sigmoid_psi = DerSigmoidCoefficient(&psi);
+
+   LinearForm int_sigmoid_psi(psi.FESpace());
+   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+   LinearForm int_der_sigmoid_psi(psi.FESpace());
+   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
+                                              der_sigmoid_psi));
+   bool done = false;
+   for (int k=0; k<max_its; k++) // Newton iteration
+   {
+      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
+      const double f = int_sigmoid_psi.Sum() - target_volume;
+
+      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
+      const double df = int_der_sigmoid_psi.Sum();
+
+      const double dc = -f/df;
+      psi += dc;
+      if (abs(dc) < tol) { done = true; break; }
+   }
+   if (!done)
+   {
+      mfem_warning("Projection reached maximum iteration without converging. Result may not be accurate.");
+   }
+   int_sigmoid_psi.Assemble();
+   return int_sigmoid_psi.Sum();
+}
+
+#ifdef MFEM_USE_MPI
+double ParVolumeProjection(ParGridFunction &psi, const double target_volume,
+                           const double tol=1e-12, const int max_its=10)
+{
+   auto sigmoid_psi = SigmoidCoefficient(&psi);
+   auto der_sigmoid_psi = DerSigmoidCoefficient(&psi);
+   auto comm = psi.ParFESpace()->GetComm();
+
+   ParLinearForm int_sigmoid_psi(psi.ParFESpace());
+   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+   ParLinearForm int_der_sigmoid_psi(psi.ParFESpace());
+   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
+                                              der_sigmoid_psi));
+   bool done = false;
+   for (int k=0; k<max_its; k++) // Newton iteration
+   {
+      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
+      const double myf = int_sigmoid_psi.Sum();
+      double f;
+      MPI_Allreduce(&myf, &f, 1, MPI_DOUBLE, MPI_SUM, comm);
+      f -= target_volume;
+
+      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
+      const double mydf = int_der_sigmoid_psi.Sum();
+      double df;
+      MPI_Allreduce(&mydf, &df, 1, MPI_DOUBLE, MPI_SUM, comm);
+
+      const double dc = -f/df;
+      psi += dc;
+      if (abs(dc) < tol) { done = true; break; }
+   }
+   int rank;
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+   if ((rank == 0) & (!done))
+   {
+      mfem_warning("Projection reached maximum iteration without converging. Result may not be accurate.");
+   }
+   const double myf = int_sigmoid_psi.Sum();
+   double f;
+   MPI_Allreduce(&myf, &f, 1, MPI_DOUBLE, MPI_SUM, comm);
+   return f;
+}
+#endif // end of MFEM_USE_MPI for ParProjit
+
+Array<int> getOffsets(Array<FiniteElementSpace*> &spaces)
+{
+   Array<int> offsets(spaces.Size() + 1);
+   offsets[0] = 0;
+   for (int i=0; i<spaces.Size(); i++)
+   {
+      offsets[i + 1] = spaces[i]->GetVSize();
+   }
+   offsets.PartialSum();
+   return offsets;
+}
+
+class BlockLinearSystem
+{
+public:
+   bool own_blocks = false;
+   BlockLinearSystem(Array<int> &offsets,
+                     Array<FiniteElementSpace*> &array_of_spaces,
+                     Array2D<int> &array_of_ess_bdr)
+      : spaces(array_of_spaces), ess_bdr(array_of_ess_bdr),
+        numSpaces(array_of_spaces.Size()),
+        A(offsets), b(offsets), prec(offsets)
+   {
+      A_forms.SetSize(numSpaces, numSpaces);
+      A_forms = nullptr;
+
+      b_forms.SetSize(numSpaces);
+      for (int i=0; i<numSpaces; i++) { b_forms[i] = new LinearForm(spaces[i], b.GetBlock(i).GetData()); }
+
+      prec.owns_blocks = true;
+   }
+
+   void SetBlockMatrix(int i, int j, Matrix *mat)
+   {
+
+      if (i == j)
+      {
+         auto bilf = static_cast<BilinearForm*>(mat);
+         if (!bilf) { mfem_error("Cannot convert provided Matrix to BilinearForm"); }
+         SetDiagBlockMatrix(i, bilf);
+         return;
+      }
+      auto bilf = static_cast<MixedBilinearForm*>(mat);
+      if (!bilf) { mfem_error("Cannot convert provided Matrix to MixedBilinearForm"); }
+      if (bilf->TestFESpace() != spaces[i])
+      {
+         mfem_error("The provided BilinearForm's test space does not match with the provided array of spaces. Check the initialization and block index");
+      }
+      if (bilf->TrialFESpace() != spaces[j])
+      {
+         mfem_error("The provided BilinearForm's trial space does not match with the provided array of spaces. Check the initialization and block index");
+      }
+      A_forms(i, j) = bilf;
+   }
+
+   void SetDiagBlockMatrix(int i, BilinearForm *bilf)
+   {
+      if (bilf->FESpace() != spaces[i]) {mfem_error("The provided BilinearForm's space does not match with the provided array of spaces. Check the initialization and block index"); }
+      A_forms(i, i) = bilf;
+   }
+
+   inline MixedBilinearForm *GetBlock(int i, int j)
+   {
+      if (i == j) { mfem_error("For diagonal block, use GetDiagBlock(i)."); }
+      return static_cast<MixedBilinearForm*>(A_forms(i, j));
+   }
+   inline BilinearForm *GetDiagBlock(int i)
+   {
+      return static_cast<BilinearForm*>(A_forms(i, i));
+   }
+   inline LinearForm *GetLinearForm(int i)
+   {
+      return b_forms[i];
+   }
+
+   /**
+    * @brief Assemble block matrices. X should contain essential boundary condition.
+    *
+    * @param x Solution vector that contains the essential boundary condition.
+    */
+   void Assemble(BlockVector &x);
+
+   void SolveDiag(BlockVector &x, const Array<int> &ordering,
+                  const bool isSPD=false);
+
+   void GMRES(BlockVector &x);
+   void PCG(BlockVector &x);
+
+   ~BlockLinearSystem() = default;
+
+private:
+   Array<FiniteElementSpace*> &spaces;
+   Array2D<int> &ess_bdr;
+   int numSpaces;
+   BlockOperator A;
+   BlockVector b;
+   Array2D<Matrix*> A_forms;
+   Array<LinearForm*> b_forms;
+   BlockDiagonalPreconditioner prec;
+};
+
+class VectorGradientGridFunctionCoefficient : public MatrixCoefficient
+{
+public:
+   VectorGradientGridFunctionCoefficient(GridFunction *gf)
+      : GridFunc(gf), vdim(gf->VectorDim()),
+        sdim(gf->FESpace()->GetMesh()->Dimension()),
+        MatrixCoefficient(
+           gf->VectorDim(), gf->FESpace()->GetMesh()->Dimension(), false) {}
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
+      if (T.mesh->GetNE() == gf_mesh->GetNE())
+      {
+         Vector V;
+         GridFunc->GetVectorGradient(T, K);
+      }
+      else
+      {
+         mfem_error("Inconsistent mesh.");
+      }
+   }
+
+private:
+   const int vdim, sdim;
+   GridFunction *GridFunc;
+};
+class FrobeniusNormCoefficient : public Coefficient
+{
+private:
+   MatrixCoefficient * a;
+
+   mutable DenseMatrix mat_a;
+public:
+   /// Construct with the two vector coefficients.  Result is \f$ A \cdot B \f$.
+   FrobeniusNormCoefficient(MatrixCoefficient &A): a(&A), mat_a(A.GetHeight(),
+                                                                   A.GetWidth()) {};
+
+   /// Set the time for internally stored coefficients
+   void SetTime(double t) { Coefficient::SetTime(t); a->SetTime(t); }
+
+   /// Reset the first vector in the inner product
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the first vector coefficient in the inner product
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Evaluate the coefficient at @a ip.
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      a->Eval(mat_a, T, ip);
+      Vector v(mat_a.GetData(), mat_a.Height()*mat_a.Width());
+      return v.Norml2();
+   }
+};
+
+
+void BlockLinearSystem::Assemble(BlockVector &x)
+{
+   Array<int> trial_ess_bdr;
+   Array<int> test_ess_bdr;
+   b = 0.0;
+   for (int row = 0; row < numSpaces; row++)
+   {
+      b_forms[row]->Assemble();
+      test_ess_bdr.MakeRef(ess_bdr.GetRow(row), ess_bdr.NumCols());
+      if (A_forms(row, row))
+      {
+         BilinearForm* bilf = this->GetDiagBlock(row);
+         delete bilf->LoseMat();
+         bilf->SetDiagonalPolicy(mfem::Operator::DIAG_ONE);
+         bilf->Assemble();
+         bilf->EliminateEssentialBC(test_ess_bdr, x.GetBlock(row), b.GetBlock(row));
+         bilf->Finalize();
+         A.SetBlock(row, row, &(bilf->SpMat()));
+         prec.SetDiagonalBlock(row, new GSSmoother(bilf->SpMat()));
+      }
+      for (int col = 0; col < numSpaces; col++)
+      {
+         if (col == row) // if diagonal, already handled using bilinear form
+         {
+            continue;
+         }
+         trial_ess_bdr.MakeRef(ess_bdr.GetRow(col), ess_bdr.NumCols());
+         if (A_forms(row, col))
+         {
+            MixedBilinearForm* bilf = this->GetBlock(row, col);
+            delete bilf->LoseMat();
+            bilf->Assemble();
+            bilf->EliminateTrialDofs(trial_ess_bdr, x.GetBlock(col), b.GetBlock(col));
+            bilf->EliminateTestDofs(test_ess_bdr);
+            bilf->Finalize();
+            A.SetBlock(row, col, &(bilf->SpMat()));
+         }
+      }
+   }
+}
+
+
+void BlockLinearSystem::SolveDiag(BlockVector &x, const Array<int> &ordering,
+                                  const bool isSPD)
+{
+   Array<int> curr_ess_bdr;
+   for (int i:ordering)
+   {
+      curr_ess_bdr.MakeRef(ess_bdr.GetRow(i), ess_bdr.NumCols());
+      b.GetBlock(i) = 0.0;
+      b_forms[i]->Assemble();
+      BilinearForm* bilf = this->GetDiagBlock(i);
+      delete bilf->LoseMat();
+      bilf->SetDiagonalPolicy(mfem::Operator::DIAG_ONE);
+      bilf->Assemble();
+      bilf->EliminateEssentialBC(curr_ess_bdr, x.GetBlock(i), b.GetBlock(i));
+      bilf->Finalize();
+      SparseMatrix &mat = this->GetDiagBlock(i)->SpMat();
+      GSSmoother curr_prec(mat);
+      if (isSPD)
+      {
+         mfem::PCG(mat, curr_prec, b.GetBlock(i), x.GetBlock(i), 0, 2000, 1e-12, 0.0);
+      }
+      else
+      {
+         mfem::GMRES(mat, curr_prec, b.GetBlock(i), x.GetBlock(i), 0, 2000, 50, 1e-12,
+                     0.0);
+      }
+   }
+}
+
+void BlockLinearSystem::GMRES(BlockVector &x)
+{
+   mfem::GMRES(A, prec, b, x, 0, 2000, 50, 1e-12, 0.0);
+}
+
+void BlockLinearSystem::PCG(BlockVector &x)
+{
+   mfem::PCG(A, prec, b, x, 0, 2000, 1e-12, 0.0);
+}
+
+// class AndersonAccelerator
+// {
+//    typedef std::__1::function<void(Vector &)> FixedPointMap;
+//    public:
+//    AndersonAccelerator(FixedPointMap f, Vector &x, const int m_=3, const double tolerance=1e-6, const int max_iteration=100): op(f), sol(x), X(m_), Gk(m_), m(m_), tol(tolerance), max_it(max_iteration) {
+//       Gk = nullptr;
+//       k = 0;
+//       converged = false;
+//    }
+//    void Solve()
+//    {
+//       X[0] = new Vector(sol);
+//       Gk[0] =
+//       for(k=)
+//       {
+//          const int m_k = std::min(m, k);
+//          FixedPointIteration();
+//          if (m < k)
+//          {
+//             X[k] = new Vector()
+//          }
+//       }
+//    }
+//    void FixedPointIteration()
+//    {
+//    }
+
+//    private:
+//    FixedPointMap op;
+//    Vector &sol;
+//    Array<Vector*> X;
+//    Array<Vector*> Gk;
+//    const int m;
+//    const double tol;
+//    const int max_it;
+//    int k;
+//    bool converged;
+
+// };
+} // end of namespace mfem
+
+#endif // end of proximalGalerkin.hpp

examples/read_mfem_mesh.m

@@ -0,0 +1,41 @@
+function mesh = read_mfem_mesh(filename)
+fid = fopen(filename, 'r');
+str = textscan(fid, '%s', 'Delimiter', '\n');
+str = str{1};
+
+%% Find Element starting point
+for i = 1 : length(str)
+    cur_line = str{i};
+    if strcmp(cur_line, 'elements')
+        break;
+    end
+end
+%% Convert Element Data
+i = i + 1;
+cur_line = str{i};
+nrElem = str2double(cur_line);
+i0 = i;
+elements = str2num(strjoin(str(i0 + 1 : i0 + nrElem), ';')); %#ok
+elements = elements(:,3:end).' + 1;
+i0 = i0 + nrElem;
+
+%% Find Vertex Starting Point
+for i = i0 + 1 : length(str)
+    cur_line = str{i};
+    if strcmp(cur_line, 'vertices')
+        break;
+    end
+end
+%% Convert Vertex Data
+i = i + 1;
+cur_line = str{i};
+nrVerts = str2double(cur_line);
+i = i + 1;
+% cur_line = str{i};
+% dim = str2double(cur_line);
+i0 = i;
+V = str2num(strjoin(str(i0 + 1 : i0 + nrVerts), ';')); %#ok
+
+%% Create SDGMesh
+V = num2cell(V, 1);
+mesh = SDGMesh(elements, V{:});