@@ -73,13 +73,18 @@ int main(int argc, char **argv) {
73
73
CeedData ceed_data ;
74
74
PetscCall ( PetscCalloc1 (1 , & ceed_data ) );
75
75
76
- OperatorApplyContext ctx_residual_ut , ctx_initial_u0 , ctx_initial_p0 ;
76
+ OperatorApplyContext ctx_residual_ut , ctx_initial_u0 , ctx_initial_p0 ,
77
+ ctx_post_Hdiv , ctx_post_H1 ;
77
78
PetscCall ( PetscCalloc1 (1 , & ctx_residual_ut ) );
78
79
PetscCall ( PetscCalloc1 (1 , & ctx_initial_u0 ) );
79
80
PetscCall ( PetscCalloc1 (1 , & ctx_initial_p0 ) );
81
+ PetscCall ( PetscCalloc1 (1 , & ctx_post_Hdiv ) );
82
+ PetscCall ( PetscCalloc1 (1 , & ctx_post_H1 ) );
80
83
ceed_data -> ctx_residual_ut = ctx_residual_ut ;
81
84
ceed_data -> ctx_initial_u0 = ctx_initial_u0 ;
82
85
ceed_data -> ctx_initial_p0 = ctx_initial_p0 ;
86
+ ceed_data -> ctx_post_Hdiv = ctx_post_Hdiv ;
87
+ ceed_data -> ctx_post_H1 = ctx_post_H1 ;
83
88
// ---------------------------------------------------------------------------
84
89
// Process command line options
85
90
// ---------------------------------------------------------------------------
@@ -106,32 +111,32 @@ int main(int argc, char **argv) {
106
111
// ---------------------------------------------------------------------------
107
112
// Create DM
108
113
// ---------------------------------------------------------------------------
109
- DM dm , dm_u0 , dm_p0 ;
114
+ DM dm , dm_u0 , dm_p0 , dm_H1 ;
110
115
PetscCall ( CreateDM (comm , vec_type , & dm ) );
111
116
PetscCall ( CreateDM (comm , vec_type , & dm_u0 ) );
112
117
PetscCall ( CreateDM (comm , vec_type , & dm_p0 ) );
118
+ PetscCall ( CreateDM (comm , vec_type , & dm_H1 ) );
113
119
// TODO: add mesh option
114
120
// perturb to have smooth random mesh
115
121
//PetscCall( PerturbVerticesSmooth(dm) );
116
122
117
123
// ---------------------------------------------------------------------------
118
124
// Setup FE
119
125
// ---------------------------------------------------------------------------
120
- SetupFE (comm , dm , dm_u0 , dm_p0 );
121
-
126
+ SetupFEHdiv (comm , dm , dm_u0 , dm_p0 );
127
+ SetupFEH1 ( problem_data , app_ctx , dm_H1 );
122
128
// ---------------------------------------------------------------------------
123
129
// Create local Force vector
124
130
// ---------------------------------------------------------------------------
125
131
Vec U ; // U=[p,u], U0=u0
126
132
PetscCall ( DMCreateGlobalVector (dm , & U ) );
127
- PetscCall ( VecZeroEntries (U ) );
128
133
129
134
// ---------------------------------------------------------------------------
130
135
// Setup libCEED
131
136
// ---------------------------------------------------------------------------
132
137
// -- Set up libCEED objects
133
- PetscCall ( SetupLibceed (dm , dm_u0 , dm_p0 , ceed , app_ctx , problem_data ,
134
- ceed_data ) );
138
+ PetscCall ( SetupLibceed (dm , dm_u0 , dm_p0 , dm_H1 , ceed , app_ctx ,
139
+ problem_data , ceed_data ) );
135
140
//CeedVectorView(force_ceed, "%12.8f", stdout);
136
141
//PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
137
142
// bc_pressure) );
@@ -192,11 +197,31 @@ int main(int argc, char **argv) {
192
197
// ---------------------------------------------------------------------------
193
198
// Save solution (paraview)
194
199
// ---------------------------------------------------------------------------
195
- PetscViewer viewer ;
196
- PetscCall ( PetscViewerVTKOpen (comm ,"solution.vtu" ,FILE_MODE_WRITE ,& viewer ) );
197
- PetscCall ( VecView (U , viewer ) );
198
- PetscCall ( PetscViewerDestroy (& viewer ) );
200
+ Vec U_H1 ;
201
+ PetscCall ( DMCreateGlobalVector (dm_H1 , & U_H1 ) );
202
+ PetscCall ( VecZeroEntries (U_H1 ) );
203
+ if (app_ctx -> view_solution ) {
204
+ PetscViewer viewer_p ;
205
+ PetscCall ( PetscViewerVTKOpen (comm ,"solution_p.vtu" ,FILE_MODE_WRITE ,
206
+ & viewer_p ) );
207
+ PetscCall ( VecView (U , viewer_p ) );
208
+ PetscCall ( PetscViewerDestroy (& viewer_p ) );
209
+
210
+ SetupProjectVelocityCtx_Hdiv (comm , dm , ceed , ceed_data ,
211
+ ceed_data -> ctx_post_Hdiv );
212
+ SetupProjectVelocityCtx_H1 (comm , dm_H1 , ceed , ceed_data ,
213
+ ceed_data -> ctx_post_H1 );
199
214
215
+ ProjectVelocity (ceed_data , U , vec_type , & U_H1 ,
216
+ ceed_data -> ctx_post_Hdiv ,
217
+ ceed_data -> ctx_post_H1 );
218
+
219
+ PetscViewer viewer_u ;
220
+ PetscCall ( PetscViewerVTKOpen (comm ,"solution_u.vtu" ,FILE_MODE_WRITE ,
221
+ & viewer_u ) );
222
+ PetscCall ( VecView (U_H1 , viewer_u ) );
223
+ PetscCall ( PetscViewerDestroy (& viewer_u ) );
224
+ }
200
225
// ---------------------------------------------------------------------------
201
226
// Free objects
202
227
// ---------------------------------------------------------------------------
@@ -205,7 +230,12 @@ int main(int argc, char **argv) {
205
230
PetscCall ( DMDestroy (& dm ) );
206
231
PetscCall ( DMDestroy (& dm_u0 ) );
207
232
PetscCall ( DMDestroy (& dm_p0 ) );
233
+ PetscCall ( DMDestroy (& dm_H1 ) );
208
234
PetscCall ( VecDestroy (& U ) );
235
+ PetscCall ( VecDestroy (& U_H1 ) );
236
+ PetscCall ( VecDestroy (& ceed_data -> ctx_residual_ut -> X_loc ) );
237
+ PetscCall ( VecDestroy (& ceed_data -> ctx_residual_ut -> X_t_loc ) );
238
+ PetscCall ( VecDestroy (& ceed_data -> ctx_residual_ut -> Y_loc ) );
209
239
if (problem_data -> has_ts ) {
210
240
PetscCall ( TSDestroy (& ts ) );
211
241
} else {
@@ -222,6 +252,8 @@ int main(int argc, char **argv) {
222
252
PetscCall ( PetscFree (ctx_initial_u0 ) );
223
253
PetscCall ( PetscFree (ctx_initial_p0 ) );
224
254
PetscCall ( PetscFree (ctx_residual_ut ) );
255
+ PetscCall ( PetscFree (ctx_post_H1 ) );
256
+ PetscCall ( PetscFree (ctx_post_Hdiv ) );
225
257
226
258
// Free libCEED objects
227
259
//CeedVectorDestroy(&bc_pressure);
0 commit comments